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An Adaptive Echo Canceller 


By M. M. SONDHI 
(Manuscript received November 3, 1966) 


A novel method 1s presented for echo-cancellation in long distance 
telephone connections. In contrast with conventional echo suppressors, the 
device described achieves echo-cancellation without interrupting the return 
path. A replica of the echo is synthesized and subtracted from the return 
signal. The replica ts synthesized by means of a filter which, under the 
control of a feedback loop, adapts to the transmission characteristic of the 
echo path and tracks variations of the path that may occur during a con-— 
versation. 

The adaptive control loop 1s described by a set of simultaneous, non- 
linear, first-order differential equations. It is shown that under ideal 
conditions, the echo converges to zero. Estimates of the rate of convergence 
are obtained. Effects of noise are discussed. The resulis of computer simula- 
tions of various alternative configurations of the system are described. 


I. INTRODUCTION 


In telephone connections that involve both 4-wire and 2-wire links, 
an echo is generated at the hybrid that connects a 4- to a 2-wire link. 
The situation at one such hybrid 1s illustrated schematically in Fig. 1. 
Here S, and S, are the two speech signals and £ is the echo of S, which 
is returned along with S, . In practice, # is an the average about 15 dB 
lower than S, , but in extreme cases may be only 6 dB lower. 

This echo has a disturbing influence on the conversation, which 
appears to increase with increasing round-trip delay.” If no steps were 
taken to reduce this echo, conversation would be seriously impaired 
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over satellite communication links with round-trip delays of hundreds 
of milliseconds. The devices used at the present time to combat echo 
are called echo suppressors. A number of different types of echo sup- 
pressors have been designed. They are all, basically, voice-operated 
switches (albeit ingenious and complex ones) which disconnect the 
return path or introduce a large attenuation in it whenever a decision 
mechanism indicates that the level of S, is large compared to that 
of S, + £. However, since # and S, both share the return path, the 
use of such echo suppressors introduces ‘‘chopping” or interruptions 
of S, during periods of double talking.* It has been shown that the 
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Fig. 1— Typical situations where an echo canceller could be used. 


degrading effect of chopping also increases with increasing round-trip 
delay.” The characteristics, advantages and disadvantages of such echo 
suppressors have been described in a number of papers.”’*"* 

It appears that improvements in such echo suppressors are not 
likely to solve the echo problem satisfactorily. Entirely different ap- 
proaches are called for. One such approach was an open loop device 
suggested by J. L. Flanagan and D. W. Hagelberger and implemented 
by J. de Barbeyrac.’ In this approach, E is regarded as a linearly filtered 
_ version of S,. The impulse response of this filter is measured by means 
of a transmitted test pulse, and a transversal filter is synthesized to 
approximate this impulse response. With S, as an input to the trans- 
versal filter, the output approximates /, and may, therefore, be sub- 


* In this paper, the term ‘“‘double-talking” will be used for the simultaneous 
presence, at the echo suppressor, of speech signals of the two speakers. 
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tracted from the return signal to cancel the echo. In this manner, 
effective echo cancellation (as opposed to suppression) is achieved 
without interrupting S,.J. de Barbeyrac demonstrated the feasibility 
and effectiveness of such a device on four- to two-wire junctions sim- 
ulated in the laboratory. 

An actual echo path is, however, not perfectly constant. Besides 
obvious step changes such as the connection or disconnection of exten- 
sion phones during a conversation, or transfer of calls via key tele- 
phones or PBX’s, there may be slow changes in gain and other fluctua- 
tions of the transfer function of the echo path. Also, economic con- 
sideration would force placement of echo cancellers at switching offices 
high in the hierarchical structure rather than at the hybrids where | 
the echoes are generated. (The number of echo cancellers required in 
the latter case would be many thousand times the number required 
in the former.) In such a situation, one or more carrier links intervene 
between the echo canceller and the hybrid. A large percentage of these 
(e.g., the N carrier) use compandors which are nonlinear elements 
with memory. Thus, what are available to the echo canceller are not 
S, and S. + # but the modified signals S/ and Si + H’ (see Fig. 1(b)). 
E’ is no longer a linearly filtered version of S{, although a linear filter 
with an impulse response dependent upon the power level of S{ could 
approximately transform S{ to E£’. Thus, for the open loop device to 
work in practice, it would seem necessary to intermittently adjust the 
transversal filter during a conversation. The transmission of test pulses 
required for such adjustments might prove quite intolerable to the 
customers. 

A proposal made by John L. Kelly, Jr. avoids these difficulties. The 
speech signal itself is used in place of test pulses and a control loop 
continuously adapts the transversal filter to take care of fluctuations 
in the echo path. 

In this paper, we will describe the system proposed by Kelly. We 
will describe various modifications of the system which simplify and 
improve it. Finally, we will report the results of tests of these systems 
by computer simulation, using artificially created echoes and also using 
two-track tape recordings made on an actual N-carrier link. 


II. KELLY’S PROPOSAL 


The adaptive control loop shown in Fig. 2 is the system proposed 
by Kelly, except for the introduction of the nonlinear function F’. 
This function is chosen to be an odd nondecreasing function with 
F(0) = 0. (Kelly’s proposal obtains as a special case when F(e) = e.) 
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Fig. 2 —Schematic of the echo canceller using a transversal filter. 
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y(t)=2(t) +n(t) 


The signal x(t) is the input speech signal (corresponding to S, or 
S! of Fig. 1). The signal y(é) is the return signal (corresponding to 
S. + EF or S§ + LH’ of Fig. 1) and is given by | 


y(t) = nt) + 2), 


where 2(é) is the echo of the input signal and n(¢) is a noise which is 
assumed statistically independent of z(t). The noise may include a 
second speech signal besides circuit noise. 

The N-tap transversal filter synthesizes an estimate of 2(¢) given by 


a(t) = d g.(é)a[t — (k — YT) 


where 7, is the delay of each section of the transversal filter. The 
control loop uses the error e(t) = y(t) — &(t) to continuously improve 
the estimate 4(é). | 

Ideally, the system should drive itself to the condition e(t) = n(é) 
(not necessarily e(t) = 0, for as mentioned earlier n(¢é) may contain 
a speech signal which must be left as undistorted as possible). Such 
ideal echo cancellation is possible with this system only if n(f) = 0 
and if 2(é) is exactly representable by passing x(t) through an N-tap 
_ transversal filter with constant (or slowly varying) tap gains. In the 
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next section we will exhibit a proof that under these conditions the 
system does indeed converge monotonically to this echo-less steady 
state. In the presence of noise the final state is one in which the tap 
gains fluctuate around their settings for perfect echo cancellation. 
The response of the system averaged over the noise ensemble will be 
shown to converge monotonically to this final state. 

The system is a first-order system and is stable for any arbitrary 
input. As is the case with most control systems, speed of response can 
be traded for immunity to noise. The constant multiplier K of Fig. 2 
allows adjustment of this trade-off. (Unfortunately, in the present case 
speed of response depends not only upon the feedback factor K but 
also upon the level and properties of the signal x(t). Thus, K can be 
adjusted to give a certain speed of response only for some average 
level of x(¢).) 

For immunity from noise, K should be made as small as possible; 
for fast convergence, 1t should be made as large as possible. We do not 
have any theory at present to calculate the optimum setting for K. 
This must be done by computer simulation and/or experiments with 
hardware implementations of the system. 


HI. CONVERGENCE 


The proof of convergence given in this section is very similar to 
a proof given by Kelly. The introduction of the nonlinear function 
F necessitates only minor modifications. However, as we shall show in 
Section V, a judicious choice of this nonlinearity can considerably 
simplify and improve the performance and implementation of the 
system. 

To simplify our discussion we introduce the following notation. We 
denote the output z[t — (k — 1)7.] of the kth tap of the transversal 
filter as x,(t). We will refer to N-tuples as vectors and consider them 
as column matrices. Thus, X(¢) will be a column matrix with elements 
x,(t), w(t), --- , ty(t), and G(é) the column matrix with elements 
gi(t), go(t), --* , gv(t). The signal 2(¢) of Fig. 2 then becomes 


| 


2) = DY anal 


= G’X. 


Here the superscript 7’ denotes the transpose of a matrix, and for 
brevity the dependence of G and X on ¢ is not explicitly shown. The 
echo z(t) will hkewise be represented as 
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a(t) = H’X, 


where H has elements h,, --+ , hy which are assumed fixed (or so 
slowly varying that their time derivatives may be neglected). Thus, 
we have 


y(t) = 2(t) + n(Z) 
= H’X + n(2) 
e(t) = y(t) — &(2) 
= R7?X + n(d, (1) 


where R = H — G. 

Before proceeding to the proof of convergence let us give an interesting 
heuristic justification for the circuit of Fig. 2. Consider a function 
C(e) such that C(e) = C(—e) and d’C/de’ = 0. C(e) is then a monoton- 
ically nondecreasing function of the magnitude of e. Let us minimize 
C(e) by varying the coefficients g,(t). If we choose to use the steepest 
descent method, we find the gradient of C(e) with respect to the g, 
and make the vector G change in the direction opposite to this gradient. 
Now 


grad C(e) = grad C(R’X + n(Z)) 
= —C'(R’X + n(t))X 
—F(R7X + n(b))X, 


where C’(-) = F(-) is the derivative of C with respect to its argument. 
To change G along the negative of the gradient we may set 

_ G = KF(R*X + n(d)X, (2) 
where KC is a positive constant of proportionality. 

By inspection, the matrix equation (2) is the equation that governs 
the dynamic behavior of the system of Fig. 2. Thus, the system is 
a steepest descent control system in the above sense. 

The proof of convergence follows easily from (2) in the case when 
n(t) = 0. Observe that since 

d d 


R=H —-G, ae ~a® 


premultiplying (2) with —2R* gives 
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ray, aor 
2R ae ~ ak R® 


= —2KR'XF(R’X). (3) 


By its definition / is a monotonic nondecreasing function and also 
an odd function. Thus, the right-hand side of (8) is always negative, 
hence R’R is nonincreasing. It is strictly decreasing whenever R*X =~ 0, 
i.e., whenever there is an uncancelled echo. Now R*R = 1? is the square 
of the length of the vector R = H — G. Thus, lz is nonincreasing, 
and as long as there is an uncancelled echo the length keeps decreasing, 
1.e., G keeps approaching H. To show that the echo goes to zero we 
integrate (3) between 0 and some time 7 and obtain 


2.0) — Bi) = K ; " R°XF(R'X) di. (4) 


As 12 is nonincreasing the left-hand side is bounded and S$/;(0). Thus, 
as tr — © we note that the integrand on the right must approach zero. 
However, the integrand is a monotonic nondecreasing function of the 
magnitude of R*X (which is the uncancelled echo). Thus, the echo 
power must approach zero.* 

If y(é) contains a noise n(t) besides the echo, then proceeding as 
before we find that 
2 R'R = —KR’XF(R’X + n(d). (5) 
From our previous discussion it follows that the right-hand side of 
(5) is negative, if and only if R7’X + n(t) has the same sign as R*X. 
As long as the magnitude of the uncancelled echo is large compared 
to n(t), this condition is met for a large percentage of time and conver- 
gence proceeds essentially monotonically as before. When the level 
of the uncancelled echo becomes of the same order as or lower than 
n(t) the convergence clearly cannot be monotonic. There will be intervals 
when n(é) is greater than R*X in magnitude and of opposite sign. How- 
ever, if the feedback gain constant K is small then G is a slowly varying 
function of time (typically K would be adjusted so that R has a ‘‘time 
constant”? on the order of 0.5 sec or so). In such a quasi-stationary 
case it is justified to assume that R*X is independent of n(t) provided 

* J. W. Sandberg has pointed out that, strictly speaking, this argument does not 
prove that | R7X|— 0 for certain pathological cases. Although these cases are of no 
practical concern, it is interesting that weak conditions suffice to rule these out. 


For example, it is sufficient that: (z) |X| and d/dt | X | be bounded and (77) the function 
F be such that eF'(e) and d/de (eF(e)) be bounded for all finite e. 
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n(t) is a wideband signal. Then it is not hard to show that the average 
of F(R*X + n(t)) over the noise ensemble has the same sign as the 
average of R*X provided only that noise has a symmetric distribution 
and thus, the system still converges in this average sense. 


IV. RATE OF CONVERGENCE 


While convergence can be proved by such relatively simple argu- 
ments, estimating the convergence rate is an extremely difficult prob- 
lem. The convergence rate depends upon the properties of X, upon 
the choice of #, and upon K and we do not have a solution to the prob- 
lem in the general case. However, we will now derive an estimate of 
the mean convergence rate for the noiseless case under the assumption 
that (3) can be averaged over the X ensemble (which is assumed sta- 
tionary) and the vector R assumed independent of X on the right-hand 
side. This assumption is justified if K is small, hence R slowly varying. 
Under this assumption, the expectation can be calculated for a variety 
of different functions F and random processes X. We will give the result 
when X is a zero mean Gaussian process with (2) F(a) = wz, and (22) 
F(z) = sgn (x) (here sgn (x) = 1 fora = 0 and —1 for z < 0). In 
case (2) 


L (R'R),, = —2Ko°(R"@R),,, (6) 


where o is the standard deviation of x;(¢) (assumed identical for all 
the x,(t)) and ® is the normalized NxWN correlation matrix of the 2;(¢). 
The angular bracket denotes ensemble averaging. In case (72) 


7 (R‘R),y — —2Ko J (R°®R) av. (7) 


These equations follow easily since any linear combination of Gaussian 
variables is another Gaussian variable. Equations (6) and (7) give 
upper and lower bounds to the convergence rates for the two cases, 
when we observe that 


Amik R S R’OR S DaaxR'R, 


where Amin and dA, are the smallest and largest eigenvalues of ©. 
If x(t) is white noise, then @ is an identity matrix with Amin = Amax = I. 
Thus, for case (2) 


VREIR Jiao XP (~KoMmnax!) 
= VR'R & VR'R |r-0 Exp (— Ko Amind) (8) 
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and in case (22) 
VRIR [2 — mo KoMhaxt 
< VR'R S VR'R |, — RE Kodnint. (9) 


These upper and lower bounds are close to each other only if ® is 
nearly diagonal (1.e., when a(t) is broadband). More elaborate methods 
can be used to estimate the convergence rate (e.g., perturbation meth- 
ods). However, although these methods would be extremely interesting 
from a theoretical point of view, they are not likely to yield much 
more insight into the convergence process. This is especially true in 
view of the fact that no satisfactory statistical description of a speech 
signal is available at present. For Gaussian noise (8) and (9) have been 
checked by computer simulation. 


V. CHOICE OF NONLINEARITY 


In the formulation of the echo suppressor problem discussed in Sec- 
tion III, the choice of the nonlinearity / depends upon the choice of 
the function C. This choice has a profound influence upon the behavior 
of the resulting system. One could set up the problem of determining 
the optimum F' which would provide, according to some reasonable 
criterion, the fastest convergence and the maximum Immunity from 
interfering noise. We do not know the solution to such a general op- 
timization problem. In any case, since convergence rate depends upon 
the statistics of the signal and noise, any such optimization would 
be practically impossible for signals as difficult to characterize as speech 
signals. We have, however, found that a number of improvements 
over the linear case result upon making F an infinite clipper. The use 
of an infinite clipper has also been suggested independently by B. F. 
Logan. 

One of the main drawbacks of making F'(e) = e (i.e., Kelly’s proposal) 
is the dependence of the time constant of the control system on the 
signal level. Equation (6), although approximate, nevertheless indicates 
that the time constant (at least for a wide band signal) is proportional 
to the signal power. A 20-dB change in signal level thus changes the 
time constant by a factor of 100. If, however, an infinite clipper is 
used the same change in signal level changes the time constant by a 
factor of only 10, which is a considerable improvement. 

Another important advantage of using an infinite clipper is the 
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considerable simplification of the circuitry. The multipliers 17,, --- Miy 
can all be replaced by switch modulators, which are far cheaper and 
simpler than broadband multipliers. 

There is another advantage in using an infinite clipper which may 
be described as follows. Suppose the system is near equilibrium and the 
echo has been reduced to a very low value. If now there is a sudden 
burst of noise (say a spurt of double-talking) then in the linear system 
the rate at which the system moves away from equilibrium is propor- 
tional to the level of the noise, whereas it is more or less independent 
of noise level if an infinite clipper is used. (In Section VI we will give 
a detailed example of this effect.) We will argue in Section VII that 
during intervals of double-talking the control loop be opened and the 
vector G frozen at its last value. However, any decision mechanism 
that would make this possible would require a finite time to make 
the decision. It is therefore important that the system should not depart 
from equilibrium too rapidly upon the introduction of a large noise 
in the return signal. 


VI. OTHER MODIFICATIONS 


It may be noted that although we started out by taking the com- 
ponents of the vector X(t) to be delayed versions of the input x(2), 
this fact was nowhere of any importance in the proof of convergence. 
All that is required is that the vector X(¢) be derived from x(t) in such 
a way that all transformations of x(t), that may possibly be produced 
by the round trip transmission path, should be representable as H’X 
with a suitable choice of H. This immediately gives us the possibility 
of generalizing the circuit of Fig. 2 to that of Fig. 3. In Fig. 3 the w;(t) 
are a set of impulse responses such that linear combinations of them 
are good approximations to most practical echo path impulse responses. 

Now there is an infinite variety of sequences of functions that are 
complete on the semi-infinite interval. One such set is the set of La- 
guerre functions which is of particular interest because it can be syn- 
thesized as a simple tapped RC ladder network. The impulse response 
of the nth Laguerre network is given by 


a TL. waa, 
L(y =e Ti (Fe ) (10) 


with the corresponding transfer function 


HO = $5 (4) a 
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_ x(t). 





e(t) =y (t)- Z(t) = y(t) =2(t)+n (t) 


Fig. 3— Generalization of the system of Fig. 2 using an orthonormal set of impulse 
responses. 


As telephone speech is limited to about 3 kHz the choice of a ap- 
propriate in the present case 1s approximately 2r xX 2000 radians 
per second, although it is not critical. 

Tests by computer simulation, to be described in the next section, 
indicate that Laguerre networks are at least as satisfactory as a tapped 
delay line for the simulation of the echo path. However, a cascade 
of RC sections would be much cheaper than a delay line. The properties 
and synthesis of Laguerre networks is described in the literature.° 


VII. COMPUTER SIMULATION 


For a computer simulation of the system described by (2), 1t was 
converted to a difference equation. Thus, if a subscript n on a quantity 
is used to denote its value at the nth sampling instant, then the equation 
simulated on the computer is 


Gret — G,, + KF(RiXn + Nn) &n- 


In one class of simulations we used filtered Gaussian noise as the 
signal and computer-simulated echo paths. Equations (8) and (9) 
appear to be very good approximations if the time constant of the 
convergence (which we may define as the time taken for 17 = R*R 
to become, say 30 percent of its initial value) is large compared. to the 
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reciprocal of the bandwidth of the input signal. If uncorrelated white 
noise is added to the echo before the system has converged then no 
convergence takes place if the noise level is about 15 dB above the echo. 
If the same noise 1s introduced after the system has converged, however, 
the balance is only slightly disturbed. A typical example will illustrate 
the orders of magnitude of these effects. The nonlinearity F was chosen 
to be an infinite clipper, and the level of the signal (which was a white 
gaussian noise) and the constant K were such that in the absence of 
noise R*R converged to about 55 dB below its initial value in 0.7 sec. 
The following two tests were performed: 


(z) The same input signal, initial conditions, etc. were used as be- 
fore, but an uncorrelated noise was added to the return signal 
at a level 18 dB above the echo. 

(22) Same as (2) except the noise was added after the system had 
converged for 0.6 sec so that 17 (hence, the echo power) was 
about 23 dB below its initial value. 


In case (7) no convergence took place and 17 hovered around its 
initial value. In case (72), after the onset of noise, lf increased slowly 
by about 3 dB in 1.5 seconds. 

For comparison the same simulations were repeated with F replaced 
by a linear function and the constant K adjusted to give a time constant 
of the exponential decay of about 0.3 sec. The noiseless case and test (7) 
gave about the same results, except, of course, that the decay was 
exponential. In test (zz) the noise was introduced after 1.1 sec instead 
of 0.6, to allow lz to converge to about 30 dB. However, in this case, 
after the introduction of the noise lz rose by about 20 dB within 0.2 
seconds. Thereafter it stayed at about this level. 

The same kind of behavior was obtained when speech signals were 
used both as input and as interfering noise and when the echos were 
generated by the computer. Of course, it is much more difficult to 
estimate time constants of the convergence when speech signals are used. 

We also used as inputs, digitized tape recordings of sentences spoken 
over an N2 carrier system. Two-track tapes were made of the input 
signal and the echo. Double-talking situations were also recorded and 
tested. For tests with these tapes, 40 to 50 delay line taps spaced 0.1 ms 
apart were used. 

With very strong echos (0-dB return loss) the system provided a 
reduction of about 20 to 25 dB as measured on a VU meter. This 
reduction took place in 0.5 to 5 seconds depending on signal level 
as discussed in Section IV. The larger variation of convergence rate 
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with signal level when the clipper is removed, was apparent in these 
tests also. | 

When a typical hybrid and a return loss of 6 dB were used the echo 
was reduced to the point of being unintelligible and almost inaudible 
even under quiet conditions. As in the case of tests with noise signals, 
the introduction of double-talking at a very high level after convergence 
had taken place produced little change in the balance. 

The fact that in the case of these recordings over the N-carrier sys- 
tem the echo could not be reduced by more than about 20 dB or so 
is undoubtedly due to the compandors used in the N-system. The echo 
canceller can provide only a linear approximation to the transmission 
path of the echo, which in the case of the N-system has nonlinearities. 

We have also simulated the echo canceller using the Laguerre expan- 
sion. In this case, the digital equivalent of the Laguerre impulse re- 
sponses were used. In terms of the delay parameter, z" = exp (—jw7’,) 
where 7’, is the sampling interval, these are 


Sh we n 
L,,(@~") ma ns (z= 4) ? 


1—az \l — az 


where a = (2 — al,)/(2 + alT’,), a being as defined in Section V. 
As mentioned earlier a is chosen so that the cut off frequency of the 
Laguerre function is about 2 or 2.3 kHz. 


VIII. DISCUSSION AND CONCLUSIONS 


We have described a new method of cancelling echos in telephone 
connections. From our theoretical discussion and simulations it appears 
that the method is feasible and can yield echo cancellation of about 
20 dB or so with a convergence time of about 0.2 to 0.5 second for 
average speech levels. This convergence time increases to 10 times 
its value for a 20-dB decrease in signal level. Convergence much faster 
than this is not possible, as then the system becomes too sensitive to 
noise and behaves erratically with normal noise level to be expected 
on a telephone connection. 

We have shown that the system would not appreciably depart from 
equilibrium even upon the incidence of double-talking. However, it 
needs, initially, a period of time in which only the echo is present in 
the return signal (of course low-level noise may be present also, but 
there should be no double-talking in this period). This initial interval 
can be as small as 0.5 seconds if the input signal is loud, but would have 
to be proportionately longer for weaker input signals. 
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It. would be advisable to break the control loop during bursts of 
loud double-talking. This need not be done by very sophisticated 
means. The following simple method should be satisfactory. Assume 
that the maximum level of an echo is 6 dB below the input signal. 
Clearly if the return signal is much larger than this it indicates double- 
talking. Rectification and integration with a time constant of about 
0.5 second gives a reasonably good estimate of the levels of the input 
and return signals. A switch could then be adjusted to open the feed- 
back path in Fig. 3 immediately following the nonlinearity F, whenever 
the input level is less than, say, 3 dB above that of the return signal. 
It is important to note that this merely prevents the gain setting 
G from changing. It does not interrupt the return path. 

We have tested our system only on an N2 carrier system (besides 
on artificially generated echos). This is a double-sideband modulation 
system in which compandors are used. There are also in use single- 
sideband carrier systems, in which no compandors are used but in 
which carrier frequency variation (during the round-trip transmission 
time) would introduce a time-varying nonlinearity. The degree to which 
this type of variation exists and its effect require further investigation. 

The delay between the input signal and its echo must be compensated 
for. This delay may be as large as 60 ms. The problem of automatically 
determining this delay and compensating for it is a challenging problem 
and is being investigated by a systems group at Bell Telephone Lab- 
oratories. They are also collecting a large sample of impulse responses 
of various connections. This information would be very useful in the 
final design of the system. For example, this will enable us to decide 
upon the optimum number of taps. Also, a fixed weighting of the gain 
vector G depending upon the statistical distribution of the impulse 
responses would improve the average performance of the system. 

The ultimate test of the system’s performance and usefulness is its 
actual use during normal long distance telephone conversations. For 
this, actual hardware must be built. Two, rather different, instrumenta- 
tions have been recently completed, one by A. J. Presti’ and the other 
by F. K. Becker and H. R. Rudin,* and it should be soon possible to 
carry out such tests. 
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T emperature Dependence of Inversion— 
Layer Frequency Response in Silicon 


By A. GOQETZBERGER and E. H. NICOLLIAN 
(Manuscript received October 28, 1966) 


Conductance-voltage and capacitance-voltage curves of metal-oxide semi- 
conductor (MOS) capacitors on n-type silicon were investigated in the 
temperature range between room temperature and 200°C. Plots of theanversion- 
layer conductance versus reciprocal temperature show a sequence of two 
activation energies: one corresponding to the temperature dependence of the 
intrinsic carrier density n; , the other to that of nz; . The low-temperature 
range is characterized by recombination-generation in the space-charge 
region, the high-temperature range by diffusion current from the bulk. The 
technique permits measurement of bulk lifetime for the two regimes and 
determination of room temperature cutoff frequency for the channel. 


I. INTRODUCTION 


Theoretical calculations of metal-oxide semiconductor (MOS) capac- 
itance show a total capacitance approaching oxide capacitance in strong 
accumulation and strong inversion.” Experimentally, it has been found 
that response time of the inversion layer can be very long.” The re- 
sponse time can be drastically shortened, however, by lateral ac current 
flow in an extended inversion layer.’'* The lateral current flow mode 
requires equilibrium surface inversion beyond the metal contact. This 
condition is usually found in p-type silicon because of the preponderance 
of positive surface charge in thermally oxidized silicon. Channel cutoff 
frequencies are then typically in the MHz range. 

In n-type silicon, charge in the inversion layer can communicate 
with the bulk under steady-state conditions only by means of genera- 
tion-recombination processes.” Inversion-layer cutoff frequencies in 
n-type silicon are normally below 100 Hz, sometimes below 1 Hz. 
These low frequencies make it difficult to measure cutoff frequencies 
and to determine the mechanism of generation of minority carriers. 
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In this study, measurements with n-type silicon were carried out at 
elevated temperature where generation is more rapid. It is thus possible 
to study the generation mechanisms and confirm the theory for cal- 
culating response time. This theory was derived by Hofstein and 
Warfield.” They consider three different generation mechanisms for 
minority carriers. These are: bulk diffusion current, space-charge gen- 
eration, and surface-state generation. Fig. 1 shows a simplified equiv- 
alent circuit proposed by Hofstein and Warfield for strong inversion. 
The inversion capacitance is fed by three parallel conductances cor- 
responding to the three generation mechanisms. Because inversion 
capacitance is large compared to oxide capacitance with which it is 
in series, it can be neglected as done in Fig. 1. 

The conductances are given for n-type bulk material by the following 
relations.” 

For surface-state response 


Gy = gBN .N per’ 'o,0, ’ (1) 


where g = electronic charge incoulombs, NV, = surface-state density/cm’, 
Ny = donor density in the bulk in em™’, 6 = q/kT, o, = capture 
cross section for holes in cm’, v, = average thermal velocity of holes 
in cm/sec, and ¥, = surface potential in volts. Relation (1) was orig- 
inally derived for a single level close to midgap. Because only levels in 
this range contribute to recombination, it is also valid for a continuum 
of surface states as is generally encountered in oxidized surfaces. 


Co 


Cp Gd 


Fig. 1—Equivalent circuit of MIS capacitor in strong inversion proposed in 
Ref. 2. C, is the oxide layer capacitance per cm2, Cp is the depletion layer capacitance 
per cm’, G,,, is the conductance arising from generation-recombination through 
surface states, mhos/cm?, G,,p is the conductance arising from generation-recombina- 
tion through states in the silicon space-charge region, mhos/cm?, and Gz is the 
conductance due to the diffusion of minority carriers from the quasi neutral region in 
the silicon to its surface, mhos/cm?. 
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Space-charge generation response: 


Nl; d 
CoS 

| ToWs 
where n; = intrinsic carrier density in em™’, 7, = bulk lifetime in 


seconds, and d = space-charge layer thickness in cm. 
Diffusion response 
2 
G, = Pee, 3 
d L,Np ( ) 
where », = hole mobility in cm’/volt-sec, and L, = diffusion length 
for holes in cm. We have further 


L, — (7 obty/B)?. (4) 


By measuring temperature dependence of the inversion-layer response, 
it is possible to determine which mechanism is dominant. Surface-state 
generation should go with an activation energy of y,. It has to be 
considered here that y, is itself a function of temperature. Space-charge 
generation has the activation energy of n; , and diffusion current that 
of n*. In the present investigation, surface-state density was made 
very small, so that only G,,p and G,z had to be considered. This was 
also done because surface-state density can reach high values close 
to the band edges.*’® This, in turn, causes considerable uncertainty 
of the value of surface potential. In the absence of surface-state effects, 
the experiments reported here showed that at low temperature space- 
charge generation dominates while at higher temperature diffusion 
current takes over. 


Il. EXPERIMENTAL TECHNIQUE 


Samples used for the measurements consisted of expitaxial layers 
of 1.5 X 10°° cm™ doping, 10 » thick, on low-resistivity substrates of 
[100] orientation. Use of epitaxial samples was advantageous because 
the measurements were not affected by series resistance in the substrate. 
Because epitaxial layers are not as perfect as regular crystals, rather 
low lifetimes were encountered. Samples were thermally oxidized in 
steam to a thickness of 1000 A. The previously described bias oxidation 
technique’ was used. In order to reduce surface-state density, the 
samples were subjected to a 30-minute annealing treatment in N, at 
350°C after an aluminum film had been evaporated.’** After annealing, 
circular areas of 3.75 X 107” cm diameter were etched out for MOS meas- 
urements. Capacitance and conductance were measured versus voltage 
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at 100 KHz and 6 KHz at various temperatures. For this purpose, 
the entire wafer was placed on a heated stage and contact was made 
to one capacitor with a wire probe. Temperature was controlled to 
+2°C. Next, depletion-layer capacitance and inversion-layer conduct- 
ance were extracted from the raw data by correcting for oxide capac- 
itance as described in Ref. 5 and 3. 


Ill. RESULTS 


A family of capacitance versus voltage curves and conductance 
versus voltage curves at 6 KHz are shown in Figs. 2 and 3. Figs. 4 and 
5 contain 100-KHz curves for the same sample. It is seen that both 
capacitance and conductance saturate in the inversion range at negative 
voltage. Due to the influence of the residual surface-state density 
small bumps appear in the depletion region. In Figs. 6 and 7, Arrhenius 
plots of the computed inversion conductance G, are presented. These 
curves were obtained from the conductance curves of Figs. 3 and 5 
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Fig. 2— Capacitance vs field plate bias measured at 6 kHz with temperature in °C 
as parameter. Sample is n-type silicon oriented in the [100] direction. Field plate 
diameter is 370 », donor density is 1.17 X 10!6 em73, and oxide layer capacitance is 
2.84 & 1077 farads/cm?. 
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Fig. 3 Equivalent parallel conductance vs field plate bias measured at 6 kHz 
with temperature in °C as parameter. Sample is the same as in Fig. 2. Conductance 
is peaked at 120°C. Dotted lines are on high and solid lines on low-temperature side 
of peak. 


at high negative voltage. The fact that both plots agree within the 
accuracy of the measurement indicates that the equivalent circuit of 
Fig. 1 is valid. The values of the activation energies also prove that 
in the surface studied here there is no noticeable influence from surface 
states. Fig. 8 contains room temperature capacitance-voltage curves 
at various frequencies. 


IV. DISCUSSION 


Hofstein and Warfield” showed that the dominant effect is most 
likely space-charge recombination (2). Surface recombination may also 
be important at relatively high surface-state densities. Because the 
sample investigated here contained very few surface states, it can be 
expected that space-charge recombination dominates. From Fig. 6 
it is seen that this is the case up to temperatures around 140°C. In 
this range, the activation energy is 0.56 eV for Fig. 6, curve (a), and 
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Fig. 4— Capacitance vs field plate bias measured at 100 kHz with temperature in 
°C as per araee: Sample is same as in Fig. 2. 


0.620 eV for Fig. 7, curve (a). The expected activation energy” for 
n; is 0.605 eV. Equation (2) can now be used to calculate bulk lifetime, 
r,, under certain simplifying assumptions given in Ref. 2. We obtain 

= 4.19 X 107° seconds. This rather low lifetime is explained by the 
fact that it refers to an epitaxial layer. 

Above 140°C a new process dominates as shown by the break in 
the 1/7 curves. This process could be either surface-state generation 
or diffusion current from the neutral part of the bulk. It can be shown 
that surface-state generation is very unlikely in this case. Surface-state 
density as determined by the conductance technique’ is varying be- 
tween 6.1 X 10’° and 3.3 X 10” states per cm’ and eV. This density 
would, according to (1), give a conductance orders of magnitude lower 
than the measured G; . It is also expected that the activation energy 
of surface-state processes should decrease because surface potential 
at constant voltage decreases considerably with increasing temperature. 

If the high-temperature points in Figs. 6 and 7 are connected by a 
straight line, they give an activation energy of 0.908 eV for 6 KHz 
and 0.935 eV for 100 KHz. This energy is lower than the expected 
energy of 1.21 eV. The discrepancy can be resolved by correcting the 
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high temperature points by subtracting the influence of space-charge 
generation as indicated in Fig. 6, curve (c). If this is done, the high- 
temperature activation energy in Fig. 6, curve (c), is 1.17 eV which 
is very close to the expected value. Fig. 7 did not contain sufficient 
experimental points to carry out the correction. 

Using (8) and (4), the high-temperature lifetime and diffusion length 
can be calculated. We find, L, = 20.1 » and7, = 1.8 X 10" seconds. 
Because the calculated diffusion length is of the order of the epitaxial 
layer thickness, it is possible that the actual diffusion length might 
be longer. In calculating the above values, a temperature dependence 
of the mobility », of JT? was used as is necessary for highly-doped 
samples. 
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Fig. 5 — Equivalent parallel conductance vs field plate bias measured at 100 kHz 
with temperature in °C as parameter. Sample is same as in Fig. 2. Conductance is 
some at 180°C. Dotted lines are on high and solid lines are on low-temperature side 
of peak, 
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The two lifetimes calculated from the two temperature ranges are 
actually expected to be equal. The linearity of the plots in curve (b) 
of Iigs. 6 and 7 indicates that there is no great temperature dependence 
of r,. A possible explanation for the discrepancy of lifetimes is that 
they are measured in different parts of the crystal. Space-charge re- 
combination occurs within 0.5 uw from the surface, while diffusion lfe- 
time is determined in the entire epitaxial layer. It is likely that a thin 
surface layer contains a higher concentration of recombination centers. 

An alternative explanation is that electron and hole lifetime are 


Gy IN MHOS 





+ (x 1073 °K7") 


Fig. 6— Equivalent parallel conductance measured at 6 kHz and a bias of —15 
volts as a function of reciprocal degrees Kelvin. The experimental points indicated 
by the circles were obtained from Fig. 3. Multiple circles at a given temperature re- 
present several runs. The solid lines are the best fit to the experimental points. 
Curve (c) is obtained by subtracting the values of G; in curve (b) from the extrap- 
olation of curve (a) at each temperature. 
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Fig. 7— Equivalent parallel conductance measured at 100 kHz and a bias of —15 
volts as a function of reciprocal degrees Kelvin. The experimental points indicated 
by the circles were obtained from Fig. 5. Multiple circles at a given temperature 
represent several runs. The solid lines are the best fit to these points. 


significantly different. In this case, 7, = (TroTpo)* would have to be 
used in (2) andr, = 7,, in (8). Under this assumption 7,, 1s calculated 
to be 10~*” second. 

By taking inversion conductance from the curves in Fig. 6 at room 
temperature, inversion-layer time constant can be accurately calculated. 
This time constant’ is tr, = Cp/G; = 2.25 X 10°° second leading to 
a cutoff frequency of 71 Hz. Fig. 8 demonstrates that a cutoff frequency 
in this neighborhood is indeed observed. 


V. CONCLUSIONS 


By measuring inversion conductance, it could be shown that the 
equivalent circuit and theory by Hofstein and Warfield is valid. The 
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Fig. 8—Capacitance vs field plate bias measured at 27°C with frequency as 
parameter. Sample is the same as in Fig. 2. 


technique applied here permits an estimate of the room-temperature 
time constant of an inversion layer by extrapolating the high-tempera- 
ture curves. This way, even very long time constants may be estimated. 
In samples having low surface-state density, like the one described here, 
only bulk generation processes are important. The temperature range 
up to 140°C is characterized by space-charge generation, above this 
range diffusion current which has a higher activation energy becomes 
more important. 
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The Charge-Control Concept in the Form 
of Equivalent Circuits, Representing a 
Link Between the Classic Large Signal 

Diode and Transistor Models 


By DANK WART KOEHLER 
(Manuscript received November 2, 1966) 


It is shown in this paper that the charge-control concept can be conceived 
as a special form of the Linvill model for semiconductors. Instead of 
mathematical tools, charge-control models become equivalent circuits amen- 
able to ordinary network analysis techniques. In the simplest form, the 
charge-control equivalent circuit for the junction transistor is fully equiv- 
alent to the Linvill and the Beaufoy-Sparkes model. For all practical 
purposes, 1t 1s also equivalent to the Ebers-M oll model. 

The charge-control junction transistor equivalent circuit combines those 
features of the other models that are tmportant for electrical engineering 
applications. It also permits the conversion between the three basic types 
of models. Because of tts close relationship to the physical processes governing 
a device, it can readily be extended to higher-order phenomena. This 1s 
usually done by expressing a Linvill-type lumped model in terms of charge 
parameters. The charge-control equivalent circuit can be useful for modeling 
a variety of semeconductor devices. 


I. INTRODUCTION 


Three basic approaches are generally used to obtain descriptive large- 
signal models for transistors and diodes, the Ebers-Moll model,’ the 
Linvill model’ and the charge-control concept? after Beaufoy and 
Sparkes. 

The Lbers-Moll transistor model’’* is based on the idea of super- 
imposing a “normal” and an “inverse”? transistor. Semiconductor 
junctions are represented by means of diodes and capacitors, whereas 
the properties of the transistor base are represented by frequency- 
dependent current sources. The Ebers-Moll transistor model is the 
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most popular of all transistor models since it lends itself most readily 
to simple “rule-of-thumb calculations.” The current relations are 
described in the frequency domain, whereas the junction voltages are 
described as functions of current in the time domain, or, as in the 
original paper, only at de. The model simulates only the effect which 
minority carrier storage exercises on the relations among the various 
device currents, but not the effect on current-voltage relations. Since 
the diode is a one-port device, no diode model of the Ebers-Moll type 
exists that could simulate carrier storage.* 

The Linvill model’’’** is almost a direct representation of the con- 
tinuity and diffusion equations for semiconductor materials. It uses 
physical rather than circuit parameters and is superior to any other 
model when it comes to incorporating second-order physical effects 
or symbolizing new structures. 

The charge-control concept stands about halfway between 
physics and circuit considerations. It has proven in the past to be 
very useful for studying storage effects in diodes and transistors, but 
appeared to be entirely a mathematical tool. Certain equivalent circuits 
have been presented’*’’’ to illustrate charge control, but, as Linvill 
phrased it, “they have little more meaning than a symbolic model 
useful for the purposes of visualizing only.”’ 

Hamilton, Lindholm and Narud compared the three models for the 
transistor in a well-written tutorial paper.’’"° They discussed the 
approximations used in deriving each model from the same physical 
background. [See also Ref. 34] In contrast to this parallel treatment 
of the three models, the following study dwells on the interrelations 
and conversions between the various models. This is illustrated sym- 
bolically in Fig. 1. 

We may call the Linvill model a physical model, the Beaufoy-Sparkes 
charge-control model a mathematical model, and the Ebers-Moll model 
an electrical model. The link between the three models is accomplished 
through a modified approach to charge-control theory: instead of 
deriving, from device physics by means of integration, mathematical 
charge-control expressions, the charge-control concept can be treated 
entirely as an equivalent circuit tool.’’ The transistor model, for ex- 
ample, is in such a form readily comparable with, and convertible 
into the Linvill and the Ebers-Moll model, provided all of these models 
are at the same level of approximation. In its simplest form, the charge- 
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_* Diode models that simulate storage and use neither the charge control nor the 
Linvill concept are usually extensions of small-signal models towards incorporating 
certain nonlinear properties. 
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Fig. 1 — Principle of derivation of transistor models and their interrelations (heavy 
lines indicate main aspect of this paper; numbers refer to conversion equations in 
the text). 


control equivalent circuit model is fully equivalent with the standard 
form of the Beaufoy-Sparkes charge-control model. But equivalency 
is usually lost, as extensions to higher-order approximations are made 
in each model. 

In this paper, we shall review the derivation of the above-mentioned 
types of models for diodes and transistors. This will be done with the 
help of a differential transmission line model. The equivalent circuit 
type charge-control concept will then be derived for diodes and tran- 
sistors. This will be followed by a discussion of higher-order approxi- 
mations, the inclusion of drift fields, and possible applications to 
other semiconductor devices. 


II. DIODE MODELS 


2.1 Mathematical Description 


As a starting point for our discussion it is assumed that the reader 
is familiar with the continuity and transport equations, describing 
current flow and carrier density in a semiconductor material. 
Continuity equations 
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—div 7,(#) = e ae +e aa (la) 

+div j,() = @ SM 4 MO — Me, (1b) 
Transport equations 

jo(t) = euplip(t) — eD, grad p(?) (2a) 

j(t) = eun,lin() + eD, grad n(?). (2b) 


jp and 4, are the hole and electron current densities, respectively. 
p and n are the hole and electron carrier densities with py and n> being 
their equilibrium values at a given temperature. FE is the electric field 
intensity. D, and D, are the hole and electron diffusion constants, 
and uw, and yu, are the respective carrier mobilities. e = +] e | is the 
value of the electronic charge. 


2.1.1 p-n Junction 


A p-n junction is described in a first-order approximation by the 
transport equation (2). The well-justified assumption is made that 
both 7, and j, are numerically small compared with the mutually 
opposing diffusion and drift currents. With the help of the Einstein 
relations 


D, = 2 Kp (8a) 
D, =u, (3b) 


and the appropriate boundary conditions one obtains the Boltzmann 
equations that express carrier densities as functions of the applied 
junction voltage v..: : 


PrO, t) = Pro EXP E bell | (4a) 


n,(0, t) = Ngo exp Ee tet | (4b) 


Here, p,(0,t) and n,(0,t) are the carrier densities on both sides of the 
junctions; pro and 7,) are the densities for v.., = 0 or, in other words, 
at points away from the junction, previously called po and 7p in (1). 
The definitions of these notations are illustrated in Fig. 2, 
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Fig. 2— Carrier density distributions in the vicinity of a p-n junction and explan- 
ation of notations used. 


In terms of excess carrier densities, (4) transform into the following 
expressions 


Dexcess(t) ie pn(0,t) —- Pno — Pal exp is va(Df a | (5a) 


Nexcess(t) = N,(0,t) — Noo = nal exp it veel} — 1. (5b) 


Together with the reasonable approximation that the hole and electron 
currents pass through the junction unchanged,* (5) uniquely char- 
acterizes the junction. 


2.1.2 p and n Regions 


The following assumptions are implied in the analysis presented 
for a p-n diode: 


(2) The p-region is so heavily doped that the electron current can be 
neglected and appreciable carrier injection occurs only in the n-region. 
(12) The problem is reduced to one-dimensional variations along the 
x axis. | 
(iz) Drift fields are neglected. (Their inclusion will be briefly dis- 
cussed later in Section 7.3.4.) 
(iv) Space charge neutrality is assumed. 


* This is not quite true for silicon diodes at low forward currents and in the 
reverse direction where recombination in the space charge layer cannot be neglected. 
With respect to some of the diode properties, especially the current-versus-voltage 
aa the discrepancy can be accounted for by changing the exponent to 
CVezxt o . 
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With these assumptions the continuity and transport equations reduce 
to 


_ 9Jp(x,t) aa te ODn(2,t) . p(x, t) — Dro : 
oe ge ee (6) 

Op, (x,t 
in(v,t) = —eD, Pe® 7 


We shall now express (6) and (7) in terms of the excess carrier den- 
SItl€S Daezcese Which we shall denote for simplicity as 7, i.e., 


p(x,t) eo Dice O)) = Drlx,t) ~~ Pro . 
Multiplying by the cross section A we obtain 


_ 01,(z,t) _—»,._: Op (a,b) p(x,t) 
oe ee (8) 
Re) S —cAD, PEA. | (9) 


These are the two equations describing the n-region. 


2.2 Differential Diode equivalent circuits 


Iiquations (8) and (9) become transmission line equations if 7, and 
p are taken as currents and voltages, respectively. (Mathematically, 
one may think of p as an analog voltage representing carrier density.) 
lig. 3 illustrates the resulting r-g-c transmission line. 

The currents in the network branches are true currents but the 
voltages associated with the nodes are analog voltages. As a reminder, 
we have labeled the nodes with encircled ‘“‘p’s’’. The series and shunt 
elements are accordingly analog resistors, conductors and capacitors 
per unit length. 

If the diode is forward biased, the junction injects carriers into the 
n-region. They diffuse across the n-region gradually recombining until, 
at x — o, all hole current is converted into electron current. Fig. 4 
shows the carrier distribution across the n-region which is equal to 
the voltage distribution along the infinitely long r-g-c line. It can 
be derived easily from (8) and (9) that, under steady-state conditions, 
the shape of the charge distribution is proportional to 


exp (=a/1,) ’ 


where 





Ly = V Dptp - (10) 
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JUNCTION N—-REGION 








= 0 dlp: = — gdx p g= 
Seo weet ap a 
(ty = Sty (1+ pot) | ss kee 8 aaa 
¥ ~ "Ee n Pno dp = — Pdx Lp T = 520, 


Fig. 3— Analog differential transmission line representation of diode model. (The 
bars indicate dimension per unit length.) 


L, 1s called the diffusion length. 7, is the hole recombination time 
constant, or hole “‘lifetime’’. 

Since the analog voltage distribution on the capacitors of the r-g-c 
transmission line is identical with the physical charge density distribu- 
tion, and since many engineers have a much better feel for the charging 
and discharging processes of such a line than for the physical process, 
the r-g-c line representation may be quite helpful as an illustration 
of the carrier injection process. In early semiconductor work, such 
r-g-c transmission lines were frequently used.”°’"**'"**'*’ No attempt was 
made, however, to attribute the physical meaning of carrier density 
to the network nodes; the junctions were represented by so-called 
K-amplifiers. These amplifiers transform the internal voltage at x = 0 


Pp (x) 





x 


Fig, 4— Excess carrier distribution in diode n-region. 
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to the external voltage with the appropriate exponential relationship, 
while not transforming current at all. _ 

Linvill’"*’*"’ has introduced new symbols for the network elements 
which relate current to carrier density. These new notations avoid 
possible confusion between analog and physical circuit parameters, 
especially voltages, and hence enable us to combine current/carrier- 
density with current/voltage networks. Fig. 5 shows such a Linvill 
model in differential form. Again we have added bars over the letters 
as it was done with the 7, g, and é in Fig. 3 to denote their dimensions 
as being ‘units per length’. 

The symbols in the models are defined as follows: 


dt.(z,t) = —H, dz p(x,t) (lla) 
di,(z,t) = —S dx ee (11b) 
dp(x,t) = —(1/H,) dx 2,(x,h), (11¢) 
where 

H, = combinance per length = eA/r, (12a) 
S = storance per length = eA (12b) 
1/H,) = eee: per length = 1/eAD,. (12c) 

diffusance 


N- REGION 





JUNCTION 











€ 


vite kT Un(i+ prone) | 


Pno 


Fig. 5— Differential Linvill diode model. Note that, in consistency with common 
transmission line notations, the reciprocal of diffusance must be used. 
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This model can be extended to include majority carriers, drift fields 
etc. The reader is referred to the literature.”"°~*° 


2.3 Integrated Diode Models 


2.3.1 Mathematical Integration 


In order to arrive at an expression for the external diode current 
from the continuity equation (8) we can integrate this expression with 
respect to x. Choosing x = 0 and x = o as the limits of integration, 
we obtain 


di,(a,t ° 4g Splat if 
= SCD) gy = [ca G4 te + + f eAp(x,t) de. (18) 


The third integral represents the total charge in the bulk material. 
With the appropriate boundary conditions 72,(0) = 12, 12,(0) 
i,(%o) = 0, 7,(%) = 7, the well-known charge-control equation’ can 
readily be obtained as 
- WO ! gt) | 
To obtain (14) from (13) the assumption must be made that A and 
7, are constant. Note that no approximations or restrictions to specific 


charge distributions are implied in (14). (They must be made, however, 
when relating the current to the junction voltage.) 


2.3.2 Lumped Linvill Diode Model 


The crudest approximation to the distributed Linvill model of Fig. 5 
is to replace the “line’’ by just one storance and one combinance’”’ 
as shown in Fig. 6. These two elements are obtained by summing, i.e., 


p(t)=p(o,t) 











LH, (t) eALp 
© p(t) Tp 
eer 

dpo(t) ~~? 

dt 


v(t) 


aes a un (+E) 


Pno 


Fig. 6—-Lumped Linvill diode model ene single-pole approximation for 
minority carrier storage. Chosen values: Ar = Ly, p(t) = p(0,t). 
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integrating, all differential storances and combinances from x = 0 
to some value Az. The value of Az is usually chosen to equal the diffusion 
length Z, . This may seem arbitrary,'” but has no effect on the terminal 
properties of the first-order model, as long as p(t) is chosen such as 
to maintain the same amount of total charge. 

The values of the circuit elements follow from (11), (12), and (5a) as 





H. = H, Ac = HL, = lt (15a) 
S = S Ax = SL, = eAL, (15b) 
2 In (1 i BO.) | (156) 
where \ 1s an abbreviated notation, used hereafter for 
€ 
Sep ees 1 
= EF (6) 


The meaning of such lumping with respect to the carrier distribution 
is illustrated in Figs. 7 and 8. The solid lines in Fig. 7 present the 
actual carrier distribution in a switching example in which a current 
pulse is assumed. As required by the transport equation, the slope 
at x = 0 is, at any time, proportional to the current. Under steady- 
state conditions, an exponential distribution is obtained. To assume 
such exponential distributions at any instant of time (dashed lines 
in Tig. 7) represents a simplifying assumption. The corresponding 


p (x,t) p(x,t) 





o | Lp 
RESIDUAL 
CHARGE 


Fig. 7 — Illustration of the (a) charging and (b) discharging process in the neutral 
bulk material. The applied signals are assumed to be forward and reverse current 
pulses. The solid lines represent the actual shape for current pulse drive; the dashed 
lines represent exponential model approximations. 
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p(x) P(x,t) 
CLUMPED WITH M=t1 











p(t)=p(o,t) 
_ LUMPED WITH 


=1 


4 
po) - —- &-- _-LUMPED WITH 
m>t 


rs 

p(o)/m--X\_------- = 7 
| 

| 

| 

| 





EXPONENTIAL 
/ 





TOTAL CHARGE =q (t)=eAL pp (0,t) 


(a) (b) 


Fig. 8— Exponential and corresponding lumped distribution of excess minority 
carriers in the bulk material of a diode. (a) Illustration of the choice of lumping 
length. (b) Time variation for m = 1; m = 1 is generally preferred in the Linvill 
model, and is irrelevant in charge models or circuit applications of the Linvill model. 


errors are negligible in all those applications where the switching times 
are large compared with the carrier redistribution times (= diffusion 
times 7, and 7, in Fig. 12). 

In the lumped Linvill model, it is assumed that the carrier density is, 
at any instant of time, constant from x = 0 toa = L,, and that it is 
0 for all x > L,. Any information on the distribution of the charge, 
especially of the slope at x = 0, as expressed in the transport equation, 
has been lost since all series elements (diffusances) are neglected. The 
only parameter of importance left is the total number of minority 
carriers and hence, the total charge. The approximation used is therefore 
equivalent to the dashed line exponential distribution in Fig. 7. 

As mentioned above and illustrated in Fig. 8(a), the length Az over 
which p is nonzero, is most conveniently chosen to equal L,. But it 
is permissible to choose Ax ~ L, if the constant value p(x) is recognized 
to be different from p(0,t); for Ax = mL, , we must choose p(t) = p(0,t)/m 
such as to yield the same total charge 


q = eAL,p(0,!). (17) 


Fig. 8(b) shows, for i = 1, the time variation of the carrier distribu- 
tions for the lumped model (solid lines) and the exponential distribution 
(dashed lines). 

As the external voltage v(t) varies, the carrier density p(0,t) changes 
accordingly. The relation between v(t) and p(0,t) has been given above 
in (15c). We shall see below that the approximation made in the lumping 
process, as discussed above, effects only v(¢) but not the current. Little 
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or no error is made whenever, and as long as the external driving 
source impedance is large. 

Since, in this and any other lumped Linvill model, all circuit param- 
eters are functions of total charges in the various sections of the semi- 
conductor we can transform the model of Fig. 6 into a charge-controlled 
model, in which all carrier densities are replaced by charges. This will 
be shown in Section 2.4.1. 


2.4 Charge-Control Diode Model 


Without invoking any of the approximations introduced in Section 
2.3.2 and illustrated in Figs. 7 and 8, the diode is completely described 
by (14) and (5a). These two equations are the basis for the classic 
charge-control theory after Beaufoy-Sparkes as applied to diodes. 
Throughout the charge-control literature, only the current appears as 
a function of the total charge but not the voltage. If we want to relate 
the Junction voltage to the charge rather than to 7p,(0,t) as it was 
done in (5a), we must make some approximation: The simplest possible 
approximation is the assumption that p(0,t) is proportional to q(Z). 
This is, for example, satisfied if the shape f(x) of the carrier distribution 
never changes, i.e., if the carriers redistribute themselves instanta- 
neously. p(a,t) is then of the form 


p(x,t) = f(x)g(). 


The shape of f(x) does not matter as long as the integral f f(x) dz 
yields the proper proportionality constant. Examples of this are the 
exponential distribution or the lumped distribution (with any arbitrary 
value of m) in Fig. 8(a). This shows the equivalency between the 
postulate of instantaneous carrier redistribution in classic charge- 
control theory and carrier density lumping in the first order Linvill 
model. 

If we now want to establish an equivalent charge-control circuit 
we must first represent (14) by corresponding circuit symbols. This is 
done in the n-region part of Fig. 9. S is the store originally introduced 
by Beaufoy and Sparkes."’’’’ To account for directionality, we have 
added a vertical bar to the store symbol in the manner of the standard 
diode symbol. The properties of S, as defined in this paper, are: 


(2) charge stored = q(t) . 
(it) current in direction indicated by arrow —>—( /-—— 


a0 
dt 
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JUNCTION | N—- REGION 


v(t) 








kT (1+ a} 


Vite e un 


PROPERTIES OF S: @ L(tt)=—— 
oR I(s)=sQ(s)-Q(0?*) 


(@) VOLTAGE = 0 


Fig. 9— Complete first-order charge-control equivalent circuit (this circuit is a 
charge representation of Fig. 6). 


(12) voltage across store = 0. 


S is often interpreted as an infinite capacitor for which 7 = dg/dt = 
d(Cv)/dt = finite, but C— © and v— 0. 


It follows from (15c) and (17) that the junction voltage is of the form 
of) = + In [1 + Ka(O)], 


where K is a proportionality factor. If we denote the steady-state 
reverse current (flowing through the diode when v(¢) is very large 
and negative) by Js we can evaluate the constant: For v ~ —© we 
obtain 





Kk-Q=-1 
and from Fig. 9 

—Ig = O/T» ° 
Hence, 

K= : 
7 I sT> 
and thus, 
a (0 |. 
vi) = x In F + ie (18) 
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Under steady-state conditions where Q/r, = I, this equation becomes 
the well-known diode equation. 

For the possibilities of incorporating the junction capacitances see 
the discussion in Section 7.3.2. 


2.4.1 Derivation of the Charge-Control Model from the Lumped Lanvill 
Model 


It represents not merely an additional proof of equivalency but also 
a good preparation for the derivation of more complex models, if we 
show” that we can derive the charge-control model from the Linvill 
model. A somewhat related modification of the Linvill model was more 
recently proposed by Beddoes.” To this end we calculate the currents 
through the elements H, and S in Fig. 6: 


t(t) = Hp, (19a) 


dp(0,t) 0. t) 


is(t) = S (19b) 


Substitution of the values for H,, S, ye p(0,t) from (15a), (15b), 
and (17) yields 


eAL, 


Tp 


ti (t) = p(0,t) = a (20a) 


dp(0,t) dato 
ts(t) = eAL, ee ae (20b) 
This result shows that the current source in the charge-control 
model of Fig. 9 represents the current 7,, through the combinance, 
and that the store S represents the current 7s through the storance. 
To find the expression for the junction, we can express p(0,t) in 
terms of q(t) by means of (17). p,9 can again be obtained from the case, 


where V — — ©, and where p(0,t) = P(0) = — pro: 
AL ; AL, 
Tai\pace =—-Isg = Ale po) | = —* pa ; 
Pp Vo— 00 


Thus, we find 





p(0, t) a q(t) Lst, 7 qt) (21) 
Dre eAL,/ eAL, Dees 


With this we can make the transition from (15c) to (18). 
2.4.2 Hvaluation of the Charge-Control Model 


The charge-control model is completely equivalent with the lumped 
Linvill model in Fig. 6; in fact, it may be considered a circuit oriented 
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form of the Linvill model. In almost all instances’’’” where the Linvill 
model is being used for circuit applications the conversion of carrier 
density into charge must be made anyhow. The charge-control equiva- 
lent circuit in Fig. 9 uses current and voltage sources plus one lesser 
known circuit element described by the simple relations 


v(t) = 0 (22a) 
(y = WO (22b) 

or, In Laplace notation 
I(s) = sQ(s) — Q00°). (22c) 


Ordinary circuit analysis techniques can be used in working with 
the model. No restriction exists with respect to the external waveforms. 
@ appears as an additional circuit parameter with additional com- 
plexity comparable to that of an additional branch current. From a 
topological viewpoint it is a branch current. This is the price to be 
paid for inclusion of the first-order dynamic storage properties. 

Junction and n-region are clearly separated in the model. Thus, little 
difficulty should arise in adding junction capacitors (dashed in Fig. 9), 
series path resistors, and leakage resistors, provided, physical knowledge 
of such effects exist. 


2.4.3 Charge-Control Model for Short-Base Diodes 


Diodes with extremely short bases do not show the exponential 
minority carrier distribution represented in Tig. 4, but rather a prac- 
tically linear fall-off (like in a transistor base except that the collector 
is now a nonrectifying contact). With reference to Figs. 3 or 5, this 
means that the distributed ‘transmission line’ is so short that the 
effect of the series diffusances H, dominates over that of the shunt 
combinances H,. The metallic contact behaves like a short circuit 
at the end of the line. 

The analogy with the r-g-c line of Fig. 8 may help the reader visualize 
the difference between the long base and the short base diode: The 
first-order approximation for the infinitely Jong line with respect to 
currents and input voltage is the parallel connection of the shunt 
resistor and the shunt capacitor; the first-order approximation for a 
very short line is the parallel connection of the series resistor and the 
shunt capacitor. In terms of the Linvill model, the short base diode 
model is obtained by replacing H, in Fig. 6 by Hz = eAD,/L and 
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L, by L. (Note that H, increases as L becomes small.) In the charge 
control model, of Fig. 9, the term 7, which represents the recombination 
time constant for the long base diode, now becomes the diffusion time 
constant. The new value 7/ can be derived most easily from the Linvill 
model as follows: 

,_ Y Q 2pOeAW Ww: 


ee ee LY | iS 


Ga,  p(O)H,/W  pOeAD,/W 2D, 








(23) 


Apart from this numerical change, the model in Fig. 9 for the normal 
diode is equally valid for the short base diode. 


2.4.4 Piecewise Innear Charge-Control Diode Model 


For many practical purposes the logarithmic voltage source relation 
can be approximated by a switch as illustrated in Fig. 10. The switch 
opens when q becomes negative and closes when gq is able to charge 
up to g > 0. A threshold voltage V,,, is connected in series with the 
forward path. If desired, the slope of the logarithmic curve 


(24) 


v=f (q) 








(a) (lo) 


Fig. 10—Piecewise linear approximation for the semiconductor junction. (a) 
Theoretical logarithmic curve. (b) Approximated curve (the dashed lines indicate the 
completion of the diode model). 
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can be added as a resistor, where J,, is an average current, which may, 
in long hand calculations, be assumed to be 3/,,,,. The saturation 
current J, must now be represented by an external current source. 


2.4.5 Application of the Model 


The above discussion of diode models serves two purposes: First, 
they form a basic understanding for deriving transistor models. Secondly, 
the diode models can be very useful in simulating dynamic effects due 
to carrier storage in diodes. 

With the piecewise linear junction approximation of Fig. 10 applied 
to the charge-control model in Fig. 9, storage time equations can be 
derived easily using Laplace transform concepts. The model has proven 
to be very useful in the analysis of step-recovery diode circuits. In the 
piecewise linear form, it can be handled without a computer, whereas, 
for the more complex models with various parasitics added, computers 
soon become mandatory. 

Switching times for step-recovery diodes are derived in Appendix A.1 
as an example of the use of the charge-control model. The equations 
obtained have been found by many authors to agree well with actual 
measurements. The normalized storage time for recovery from an 
infinite ON-pulse according to (97) is plotted in curve a of Fig. 11 as 
a function of the reverse-to-forward current ratio according to the 
relation 


DS eln (1 4- zz), (25) 
Ip } 

When applying this result to an ordinary diode with homogeneous 
doping profile, one must be aware of the implied approximations: 
(t) The single-section approximation in the model does not affect any 
mutual relationships between currents and charges, but represents 
approximations with respect to the junction voltage. As the amount 
of stored charge is reduced considerably in the diode, the junction 
voltage decreases noticeably. (77) As the carrier density near the junction 
becomes extremely small, the voltage reverses sign and the diode 
impedance, at some point, becomes comparable with the external 
source impedance. The ideal current source assumed in (97) ceases 
to exist, and instead of the step-recovery, as given by the model, a 
long tail in the current response results. 

From either one of the two differential models in Figs. 3 or 5, we can 
calculate the time in which the carrier density at x = 0, and hence 
the junction voltage, reaches zero. Such a calculation yields the relation 
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leet 


“0.01 0.02 0.05 Of 





Fig. 11—Comparison of diode storage times as functions of the driving ratios. 
(a) Single lump model; 7 = time when charge is fully depleted. (b) Differential 
model; 7 = time when excess carrier density at x = 0 reaches zero. 


originally derived by Lax and Neustadter”’ 


T 1 
erf J = Tn (26) 
F 


l+7 


This relationship is illustrated in Fig. 11(b). Since curve (b) represents 
only the storage phase but not the very long tail of the recovery, the 
values are much smaller than those in curve (a) in which some sort 
of “effective total recovery” is represented. The difference is most 
remarkable at strong relative reverse drives where the carrier distribu- 
tions on the lines differ most from the steady-state distributions. 

If it becomes necessary to incorporate the tail of the recovery into 
a lumped diode model, the double z-extension described below may 
prove adequate for most applications. 
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2.4.6 Higher-Order Approximations 


Bearing in mind how the model originated as an approximation to 
the differential transmission line or as just another form of the lumped 
Linvill model we can now understand how higher-order approximations 
are to be obtained. 


Fig. 12 shows the example of a z-approximation for a diode. The 








WHERE ig (t)= a 2 nee = = _ gett 
_ kT q(t) 
v(t) = “@ Ln (v4 a | 


Te 
T, (4 -1)+T, 
CHARGE CONSERVATION CONSTRAINT ON a,b,C: 
cb (i-a) Lp? 


(a+b-1) Dp To _ 


Ty (FROM DC CONSIDERATIONS) = 


WHERE C IS MOST APPROPRIATELY CHOSEN TO BE c= a 


(b) 


Fig. 12— Higher-order, z-Approximation of diode charge-control model. (a) 
Charge approximation. (b) Corresponding model. 
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charge is broken up into two parts g, and q.. The diffusance between 
the two stores controls the redistribution of the charge. Such a structure 
provides a better representation of the junction at the higher frequencies 
or at higher speeds than the model of Fig. 9, since the junction voltage 
is now a function of only that part of the total charge which is close 
to the junction. The model simulates recovery tails. It also permits 
the simulation of variations in recombination time along the z-axis. 
Fig. 12 assumes two different recombination times 7, and 7.. The 
t = f(q) relation then becomes 


= d(qi + qe) na q2 
71> di + 7; + T5 (27) 


[which reduces to (14), if one assumes 7) = 79]. 

Three additional time constants 7,, 7, and 7, appear in Fig. 12. 
They depend on the choice of the sections al, and bL, over which the 
shunt elements are integrated and on the choice of the section cL, over 
which the diffusances are integrated. The three degrees of freedom 
reduce to one, however, if one considers that (2) the total charge must 
be conserved by the lumped approximation, and (zz) in a multisec- 
tional approximation the diffusances are most appropriately lumped 
over sections cl, which extend between the centers of the charge 
sections. The corresponding relations are given in Fig. 12; derivations 
have been omitted. 


Ill. LARGE-SIGNAL TRANSISTOR MODELS 


In complete analogy to the diode models, we shall now compare 
the various junction transistor models and establish the charge-control 
model in the form of equivalent circuits. The Ebers-Moll concept, 
which was found not to be applicable to dynamic diode description, 
will now enter the ‘‘competition’’. 

In order to dwell on the philosophies underlying each concept we 
shall, at first, limit ourselves to diffusion type junction transistors, 
neglecting again drift currents and secondary effects such as base-width 
modulation. All derivations will be carried out for pnp transistors; 
but, of course, everything will be correspondingly valid for npn tran- 
sistors. 


3.1 Differential Transistor Model 


The most rigorous of all the equivalent circuits describing a Junction 
transistor, as defined by (5), (8), and (9), is the differential model shown 
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E fe p(0,t) p (W,t) C 
VWs Fly} — 
E 
dH-{|_ | | |ds 
x 
ee rp 2=W 
Ve(t) =f [p(o,t)] ls NN ee ad [P (w,t)] 


Fig. 18 — Differential Linvill model for the transistor with drift fields neglected, 
pnp version shown. 


in Fig. 13. Linvill notations comparable to the diode model in Fig. 5 
were chosen. (If so wanted, the model could also be drawn with the 
notations used in Fig. 3 resulting in an r-g-c line and two K-amplifiers 
at both ends.) 

The base section of the transistor model is only a very short ‘“‘trans- 
mission” line when compared with the “infinitely long’ diode n-region 
of the normal diode. Instead of 100 percent recombination, as found 
in the diode, the transistor must have as little recombination as possible 
in order to achieve high gain. Fig. 14(a) shows a steady-state charge 
distribution under normal forward operation, and Fig. 14(b) shows 
the distribution for the case where both junctions are emitting, 1.e., 


p (o,t) 
p (0,t) 


p(x,t) 
x x 
0 W 0 W 


(a) | (b) 


Fig. 14—Excess minority carrier distribution in the transistor base under (a) 
normal and (b) saturated operation. 
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in saturation. Under normal operation the collector acts as a ‘‘charge 
short circuit”? for the line. For high-gain units, the slope is almost 
a straight line; at « = 0 it is proportional to J, and at x = W propor- 
tional to lg. 

The general case is that of Fig. 14(b) where both junctions are 
emitting and p(W) =~ 0. Any section of the base region can be de- 
scribed analogously to a four-pole using the definitions given in Fig. 15. 
Note that there are no nonlinearities in the base section. 


Te(s) = AuPils) + Ai2P2(s) (28a) 
Te(s) = AnP.(s) + Ao2P2(s). (28b) 


By using complete analogy to standard transmission line theory, it 
can be shown that with the use of (10) and (11) one obtains for a 
homogeneous section Az 




















Tz(s) = i coth y Av — oe cosech y Ax (29a) 
I¢(s) = us. cosech y Ax — a coth y Az, (29b) 
where 
ae - \z 1 ; ST (80) 
VS V1 + sr. (31) 


In the general case, the base is not homogeneous, which means that 
Z and y will vary along the line. 
The junctions are described by the time relations 


pi(0,t) = Prolexp {drv,(Z)} — 1] (32a) 
pW ,t) = prolexp {rv-(A)} — 1]. (32b) 





Fig. 15—-Symbols and polarity conventions defining the four-pole description of 
the ransintor base. 
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Fig. 16 — First-order lumped Linvill transistor model. 


IV. THE LUMPED LINVILL 7-MODEL 


In the literature, any lumped approximation to the transmission 
line model presented in Fig. 13 is referred to as a “lumped Linvill 
model” or simply ‘‘lumped model’. The most common form is the 
w-model. With respect to its current properties and steady-state voltages, 
this form will prove to be equivalent to the commonly known ‘Ebers- 
Moll model’, if both models are taken to be in the form of first-order 
approximations. 

By integrating all differential diffusances over a length Ax = W, 
all emitter sided differential combinances and storances over a length 
Ax = W,, and all collector sided differential combinances and storances 
over a length Av = W.(= W — W,) one obtains the circuit shown 
in Fig. 16. Nonsymmetry has been taken into account by using different 
recombination times 7, and 7. and different cross-sectional areas A, 
and A, on the two sides. Note that the latter represents an extension 
from the one-dimensional carrier flow and as such an example for the 
reduction of multidimensional effects to a one-dimensional model. 
Area A is some average cross section effective for the diffusion process. 
H, 1s the diffusance, the H.’s are the combinances, and the S’s are 
the two storances. 
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The four-pole equations describing the transistor base are obtained as 


11(0) = *4P2 po) 1 4 AH 1 + on) | AP PY 





= HPA9| 1 a a ( a ‘| — HP.(s) (38a) 











o(9) = 72 py) — 422 pysy [1+ Sa Ws 4 + on) | 


£)], (33b) 

The junctions are described as 
pit) = p(0,t) = prolexp {rv.()} — 1 (34a) 
pt) = p(W,t) = prolexp {rv.(t)} — 1]. (34b) 


A constraint has to be satisfied: Under equilibrium conditions, the 
total charge in the base must equal that in the two sections, 1.e., 


I 





H,P,(s) — HPA) 1 4 A ( 


Ww 
PAWEP: Ach WP? Sed | P(2) dx © 4eAWIPO) + P(W)]. (35) 


The approximation holds for high-gain units. For this case the base 
volume sections are equal, ie, A4,\W, = A.W. = Z3AW. For low-gain 
units (84) must be modified: The terms p,(t) or p.(t), or both, must 
be replaced by p,(t)/m, and p.(t)/m., respectively, whereby the m’s 
are constants >1, similar to m in Fig. 8. 

Equation (83) represents one of several possible approximations to 
(29) with the additional property of nonsymmetry being added. 

Higher-order approximations of a. lumped linear model are obtained 
by representing the base of width W by more than the two sections 
W, and W,. 


V. THE EBERS-MOLL TRANSISTOR MODEL 


The focal point of the Ebers-Moll model is the two-port description 
of the base. Such a description has been given in (28) and (29), and 
is permissible because of the linearity which exists between currents 
and carrier densities in the base. Nonlinearity exists, however, in the 
relationship between carrier densities and external voltages according 
to (32). Since linearity allows the use of the superposition principle, 
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the total current can be conceived as consisting of the superimposed 
contributions of the currents injected by the two junctions. 

When put into the form of an equivalent circuit, the Ebers-Moll 
model shows the superposition of a normal transistor (subscript NV) 
and an inverse transistor (subscript J). In Fig. 17(a) the lower diode 
and current source represent the normal transistor and the upper 
elements represent the inverse transistor. Each junction is represented 
by a diode, a fraction of the diode current is collected by the other 
electrode. The ratios of collected currents to emitted currents are called 
ay and a, for normal and inverse operation, respectively. The general 
frequency behavior of the a’s can be calculated for a homogeneous 
base from (29), (80), and (81) as 





‘ EVe (t) 





ler (t) = lero [exe(SE -') 


4B (a) 





—arIc(s) Ld tt) =Ic9| exe (ee -i) 

















OB 
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Fig. 17— The two forms of the Ebers-Moll transistor model: (a) direct representa- 
tion of the idea of superimposing a normal and an inverse transistor, (b) collecting 
current sources as functions of the electrode currents. The junction saturation currents 
in (b) are identical with the open-electrode diode saturation currents. 
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One a _ cosech y Ax _ i 
L snput() Poutput =0 coth Y Ax cosh Y Ax 
if 2D,7 
y Ax ar ~ Dir + Ax? + sr Ax” (36) 


ean 


But symmetry does not exist in a practical transistor. The constants 
in (86) are therefore, different for ay and for a;. Equation (36) can 
be rewritten under this consideration in the well-known form 


ONO 

a) — TF shaun 
aro 

ay(s) = ea ee (37b) 


The relations between the constants in (37) and the physical param- 
eters (corresponding to the constants in (36) modified for the non- 
symmetrical case) will be derived in Section 5.1. 

On account of their nonlinearity, the junction diodes must be de- 
scribed in the time domain. In their original paper, Ebers and Moll 
defined only a de relationship between voltages and currents. This 
would restrict the use of their model to piecewise linear analysis. 
But the model can be made more general®* by postulating that the 
v = f(z) relation be valid at all times, as indicated in Fig. 17. 

In either case, an important property of the semiconductor junction ‘ 
is lost: Voltages and currents appear as being directly related instead 
of being related indirectly through current density or charge. This can 
best be illustrated by anexample. If a forward current through a junction 
is suddenly replaced by a reverse current the voltage actually does not 
reverse sign until the excess carrier density at the junction is reduced 
to zero. According to the Ebers-Moll model, voltage and current always 
change polarity together. As mentioned before, it is for this reason that 
for a diode, no dynamic model of the Ebers-Moll type exists that 
would represent charge storage effects. In addition to this shortcoming, 
the feature of mixed time and frequency domain characterization is 
undersirable if the model is to be used in its nonlinear form, say on 
a computer. 

The Ebers-Moll model was originally presented in a form, shown 
in Fig. 17(b), which differs slightly from that in Fig. 17(a). Both ver- 
sions have been used throughout the literature over the past years 
and very few authors’’*° have clearly pointed out the difference between 
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them. In Fig. 17(b) the collecting currents are a times as large as the 
total emitter and collector currents, respectively. A simple calculation 
shows that the two versions are formally equivalent, if the relations 


Lur(s) = << (38a) 


and 


I¢(s) 
vi or(s) t= aty(s) ay(s) (38b) 
are satisfied. A glance at the equations for the voltage sources in Fig. 17 
reveals that the two versions could not be completely equivalent, 
unless either Jg7o and Igo or Iz, and Ig, would be considered frequency 
dependent. Due to the approximative nature of both models, this is 
normally not done. 
From both Fig. 17(a) and (b) the respective four-pole equations, 
on which the model is based, can readily be derived in terms of elec- 
trical parameters: 


Ii(s) — ar(s)IG(s) 


22.050) (39a) 


Ix(s) = Ler(s) — ar(s)Ier(s) = 
LO =a@ie@ 146) gee oe (39b) 


After substituting the expressions for the junctions one obtains for 
the steady-state case the well-known Ebers-Moll equations 


> T36 T\ _ arol co _ 
Iz = ere [exp AV.) — 1] ee [exp (AV.) — 1] (40a) 
feo __ texn (xV,) — 1]. (40b) 


pos Oil = 


] — AnoM@yro ] — Ano&ro 


5.1 Comparison Between the Ebers-Moll and the Linvill Model 


Comparing (40) with (83) and (34) for the steady-state solution 
leads to the following relations: 


— — ee — re 9 


W 1 — ayoQzo 1 — ayoetro 


eA D, Pro ayol x0 arol co (41) 


A corresponding comparison for the ac case would yield the same 
expression as in (41), except that ay» and a;) would have to be replaced 
by their frequency dependent forms. Since the left side term of (41) 
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is frequency independent, no rigorous equality exists between the 
Linvill model and any of the two versions of the Ebers-Moll model 
under ac conditions. In the Ebers-Moll model, the junction voltage 
is a function of the total diode current [being different in the two 
versions of Fig. 17(a) and (b)]; in the Linvill model it is only a function 
of the resistive component of the diode current in Fig. 17(a); this 
component equals the current through the combinance which is propor- 
tional to the carrier density p. It can be shown that the correct solution 
in which the junction voltage is a function of the carrier density directly 
at the junction, lies between these two cases but much closer to the 
lumped Linvill simulation. The discrepancy, mentioned here, affects 
only the junction voltages and does not appear in many analyses 
that use piecewise linearity. 

ay(s) and a;(s) can be expressed in terms of the physical parameters 
by comparing (89) and (83) separately for the normal operation (I¢r=0) 


and for the inverse operation (zr = 0). Subsequent conversion of 
the a’s into 6’s yields 
= ay(s) _ Byo - A7,D,/A,WW, 
Of) = Tow) 4 1 + s7, (42) 
WBN 
os ay(s) os Bro = Ar.D,/A,WW>. 
Br) = 1 — a;(s) — 1+ 8 a 1 + s72 (43) 
Wer 
where 
—_ Wan es 1 
nee oe 1 -} Byo Ty (44) 
and 
a -) eee lL (45) 
Oe 1+ Bm 7 
By definition we shall call in later sections 
T, = Ten (46) 
and 
T2 = Ter (47) 


5.2 A Better Approximation for the a Frequency Dependence in the Ebers- 
Moll Model 


Pritchard*’ has first suggested that a better approximation for the 
3-dB cut-off points of the a’s or #’s are obtained if one inserts a factor 
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1.22 into the corresponding equations, 1.e., 


1 1 . 
(SS eS (48) 
cosh y Ax 1+ 1.22% 


Weut-off measured 


This can readily be calculated from the fall-off behavior of the cosh 
expression while assuming Byp > 1. 

The same factor 1.22 appears in the corresponding expressions for 
a, , By and 6, . It is evident from (48) that this problem can be reduced 
to a matter of defining w,,. For less ideal transistors the factor is 
usually between 1 and 1.22. 

Higher-order approximations to the hyperbolic function commonly 
use two pole expressions or delay-producing excess phase terms. 


VI. THE CHARGE—-CONTROL TRANSISTOR MODEL 


In analogy to the diode charge-control model we can establish a 
charge-control equivalent circuit for the transistor. To that end, we 
want to express all parameters in terms of the charge in the base. 

Three approaches appear feasible: A lumped Linvill model can be 
labeled in such a way that all elements appear as functions of charges 
rather than integrated charge densities of the form pAzx . The two are 
proportional; the proportionality factors are of the form ‘‘electron 
charge times area’. Most of the special circuit components of the 
Linvill model become current or voltage sources in the charge-control 
version. This procedure of converting a Linvill model into a charge 
control model can readily be applied to higher-order Linvill models. 

A second approach is to use the Ebers-Moll principle of superposition 
whereby two charge-control diode models plus the corresponding col- 
lecting currents can be joined to form the transistor model. This ap- 
proach is essentially limited to’ the first order of approximation. Two 
seemingly different, but fully equivalent and easily convertible models 
result. | 

The third and classic approach to charge-control theory, originated 
by Beaufoy and Sparkes,’ is basically mathematical. Through integra- 
tion of the continuity equation the carrier density as a variable is 
replaced by the total charge in the base. Certain simplifying assump- 
tions have to be made to obtain a relation between currents and charges. 
In essence, these assumptions are equivalent to the approximations 
implied in the first-order Linvill and Ebers-Moll models as well ‘as in 
the first-order charge-control equivalent circuits to be described below. 
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Some equivalent circuits have been presented in the literature, but 
they were less rigorous than the circuits described below in the sense 
that they cannot be used as complete networks. Additional knowledge 
of the physics of the device is required to use these models. Extension 
to higher-order models in the Beaufoy-Sparkes approach is accom- 
plished through increased physical and mathematical complexity and 
not through more complex network topology as in the Linvill model 
or the charge-control model to be described. 


6.1 The x-Version (Base-Controlled Version) of the Charge-Control 
Equivalent Circuit 


In the lumped Linvill z-model of Fig. 16, the base charge distribution 
is approximated by two levels of carrier density. This is illustrated 
in Fig. 18. p,(é) is constant over the length Ar = W,, and p,(é) is 
constant over the length W,, where W, + W. = basewidth W. The 
total charge in the two sections follows with (84a) and (34b) as 


q(t) = p,(t)WyeA, = proWieA,[exp {rv,(4)} — 1] (49) 
g(t) = p.(t)W.eAs = ProW.eA.2[exp {rv,(t)} — 1]. (50) 


Using the definitions of the elements given in Fig. 16, one can calculate 
from the Linvill model in Tig. 16 the currents through H.,, H..2, 
and H, and obtains 











a(t) = = H.,p,(t) = AW a p(t) = aud = ast (51) 
ae tat a aa pO ue = ao (52) 





Fig. -18— Excess carrier distribution in the transistor base as used or implied in all 
eerie transistor models. (a) General case p: ¥ p(0), pe ¥ pt W). (b) Commonly 
used choice p; = p(0), pe = p(W). 
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t3(t) = H[p,(t) — po(t)] 
= APs p(y — vo) = Pek (BO. — HO). 30 


Using the more familiar Ebers-Moll notations and the relations found 
earlier in (42) through (47), 7; can be expressed as 


is(d) = © gu() — ™ ald (53b) 


Thus, the three current sources in the ceils model are found 
and related to the Linvill model by means of (51) through (53). 

The remaining two branch currents 2. and 2, are obtained from 
Fig. 16 as 


dpi(t dqv(t 
in(t) = eds W, SAD — Sew) (54) 
i()) = eA.W, ie = oe (55) 


These equations describe two stores Sy and §; , whose properties have 
been described in Section 2.4. 

The conversion between the two models will be summarized and 
further discussed in Section 6.4. 

The voltage sources for the junctions follow from (84), (51), and 
(52) as 


dv, (t) = In ( + ng = In (i + — wld) (56) 
dv,(t) = In (1 a ert) = In (1 4 ho). (57) 


With the help of (41) through (47) that link the constants used in 
the Ebers-Moll model to those in the Linvill Model, (56) and (57) 
can be rewritten as 


= q(t) |» ___ + Bwo | 

rv(Z) =n c T TBN - Tno/(1 7 | (58) 
uO 1 + Bro | 

Nve(i) = in [1 rol ee BI a co/(L — anoero) oe) 


(The reader may prefer to derive the constants directly from the 
steady-state Ebers-Moll model in (40) by considering the limiting 
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cases v, = 0 and v, = 0.) With the addition of (58) and (59) the equiv- 
alent circuit in Fig. 19(a) is completely defined. 

It is customary and useful in charge-control work to define additional 
parameters rzy , Ten , Ter, Tcr. We define their relationship as follows: 


Ten _ TEN _ 
Byo QNno “ae ey) 
TBI TCI 
ee re 61 
Bro QI0 sic ( 
Since, as usual, 
Byo 
= 62 
oN LF Bye a 
and 
Bro 
Qo = 1 | Bro ’ (63) 
it also follows that 
1 | 1 
a a p (64) 
TEN TBN TCN 
and 
1 1 1 
== 4. (65) 
TCI TBI TEI 


(The classic definition of these time constants after Beaufoy-Sparkes 
will be discussed in Section 6.5.) The subscripts B, H, and C stand 
for base, emitter, and collector, respectively. The subscripts N and J 
on the time constants and on the charges have been chosen to indicate 
the normal and inverse transistor operation. Many authors’’’’’’’ use 
F (forward) and R(reverse) instead of N and J. Since F and R are 
commonly reserved for diode forward and reverse currents, and since 
such currents can flow in each of the two junctions, the different nota- 
tions N and J, as proposed by Ebers and Moll, appear more appropriate. 

In Appendix B, the notations for the stored charges and the time 
constants used in this paper are related to those used in a recent book 
published by the Semiconductor Electronics Education Committee;” 
they are also compared with the notations and definitions used by 
Beaufoy and Sparkes. 

The additional time constants do not add any additional degree 
of freedom. But it is advantageous to use ‘‘base’ notations when 
controlling base current, 1.e., in common-emitter or common-collector 
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Fig. 19 Charge-control equivalent circuit for transistor in first-order approxi- 
mation, shown in two equivalent and convertible forms: (a) z-version, Linvill type, 
(b) T-version, Ebers-Moll type. 
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connection, and to use “emitter” and “collector” notations when 
emitter and collector forward currents are injected, such as in common- 
base connection. 

The model obtained in Fig. 19(a) maintains the most valuable 
property found in the Linvill model, namely the close relationship 
between the physical processes and the circuit elements. For example, 
Tay and rz, are the recombination times on the emitter and collector 
side, respectively, and rey and rz; are the diffusion time constants 
for the charges injected from the two junctions. Junctions and base 
are represented by individual sections within the equivalent circuit. 
This separation makes it easy to expand the model and to take other 
effects into account. 


6.2 The T-Version (Emitter-Collector Controlled Version) of the Charge- 
Control Equivalent Circutt 


In complete analogy with the derivation of the Ebers-Moll model 
in Fig. 17(a) we can take two diode charge-control models back to 
back and add current sources on the collector and emitter side, which 
ATe ayo and ayo times the diode currents. 

For the simulation of the junctions, we are left with two alternatives: 
One is to convert the corresponding expressions in the Ebers-Moll 
model in Fig. 17(a) into charge functions; the other is to use the expres- 
sions in the charge-control 7-model (which are equivalent to the Linvill 
model), but replace the 6-notations by a-notations according to (60) 
through (63). The first-mentioned alternative for simulating the voltage 
sources would amount to simply substituting the diodes from Fig. 17(a) 
for the voltage sources in Fig. 19(b). The property of charge control 
would not be simulated. The second procedure is therefore chosen; 
it yields 


v(t) = in (1 + ee ee 0 ces) (66) 
v(t) = xin ( = ee at (67) 


Thus, the equivalent circuit in Fig. 19(b) is obtained. 

As far as the current relations in the models are concerned, the 
main difference between the charge-control 7-model and the Ebers-Moll 
model is that the frequency dependence is simulated by a mathematical 
expression in the Ebers-Moll model, and by an additional network 
branch in the charge-control model. This is analogous to the option 
existing in small signal models where one can represent the frequency 
dependence either with an appropriate RC circuit, holding a, frequency 
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independent, or alternatively, with a frequency dependent a@ in the 
collecting current source. 

The equivalency of the charge-control model with the Ebers-Moll 
model exists only for the relations between the currents. It can be 
shown readily that the following relations must be satisfied to establish 
equivalency: 





(2) ie are : (68) 
oF 1 
(12) Tor = ae (69) 
(iit) Thy = i aa 1 + Bwo (70) 
WBN Wan 
(iv) Tar = al aa 1+ Bn. (71) 
Wr Wal 


All w’s must be replaced in these equations by the corresponding 
w/1.22 if the w’s correspond to the measured 3-dB gain fall-off points, 
and if the better approximation mentioned in Section 5.2 is to be 
included in the Ebers-Moll model, provided the particular transistor 
follows the underlying theory well enough. 


6.3 Conversion Between the Two Proposed Charge-Control Models 


The identity between the two charge-control models, presented in 
Figs. 19(a) and (b) can best be proven by converting one model into 
the other. 

To convert the z-model into the 7-model one first adds a branch 
current 2; both into and out of the base point B’ and splits 2, up into 
its two components. The resulting circuit diagram is shown in Fig. 20. 
The two parallel current sources proportional to gy(¢) on the left side 
can then be combined into one current source. Likewise, the two current 
sources proportional to q;(¢) on the right side can be combined. If 
with the help of (60) through (63), one now relabels all current sources 
in terms of ay and a; instead of 8y and 6; and extends the upper current 
sources beyond the voltage sources, one obtains the model in Fig. 19(b). 


6.4 Summary of the Conversion Equations between the Linvill and the 
Charge-Control Model 


6.4.1 Conversion Equations for the First-Order Transistor Model 


In (49) through (59), the charge-control 7-model was derived from 
the Linvill model. With the help of the defining equations for the 
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q(t) a(t) 
Pro TBI Br, Ter 





Fig. 20 —Intermediate step used in the conversion from the z to the 7 charge- 
control model, demonstrating equivalency between these two models. 


Linvill model elements, the constants in the charge-control model can 
be calculated as a function of the Linvill combinances, storances, and 
diffusances. For the relation between the Linvill 7-model of Fig. 16 
and the charge-control +-model of Fig. 19(a), such calculations yield 


S; 








rw = 7 (72) 
rn = Fe (73) 
rox = (74) 
Ter = 7 (75) 
Teo = Dan ttle + a Z H.)Hs (76) 


Note that (72) through (75) reveal that the five parameters in 
the Linvill model lead to only four parameters in the charge-control 
model. The one degree of freedom that is lost in the charge-control 
model is the conversion factor from current to carrier density; con- 
version of the charge-control model into a Linvill model is only possible, 
if one of the five Linvill parameters is known. This is tantamount to 
saying that one needs some information on the geometry of the device 
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such as the value of one, or in low-gain units, both of the two base 
volume sections A,W, and A.W.. 


6.4.2 Conversion Between the Linvill and the Charge-Control Model for 
an Arbitrary Number of Base Sections 


In higher-order approximations for diodes or transistors, the param- 
eters of the Linvill and the charge-control models, as defined in Fig. 21, 
are related by the equations 


Si S. S, 






















i Vas Ho (ha As ™ = H., (78) 
8, va 

T12a — ioe Tuva — Haus (79) 
So S, 

T12b = Ae Tyuvb = Haw (80) 

(81) 

(82) 


{ qi 
en + — 
¥, cun( i) 


Fig. 21—A Linvill and a charge-control equivalent circuit for a junction and part 
of a multisectional n-region, with indication of the notations used in converting one 


model into the other. 
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6.5 The Transistor in Saturation 


In all lumped transistor models (Linvill, Ebers-Moll, or charge- 
control type) the charge in the base is explicitly or implicitly broken 
up into the charge gy injected from the emitter under normal operation 
and the charge q; injected from the collector under inverse operation, e.g., 
in saturation. This was illustrated in Fig. 18. 

When the transistor is overdriven into saturation with a base current 
larger than I¢ s.t/Byo , the two stores gy and q; do not change by exactly 
equal amounts, 1.e., 


dqn dqr 
dip excess - dtp excess 
where, by definition, 
F ‘ I 8a 
Up excess = t3 — gE (83) 
0 


This is illustrated in Fig. 22. It can be calculated from any of the two 
models of Fig. 19 that, under steady-state conditions, the excess charges 
in the two stores are related to the excess base current by the expressions 


AnoT 
AQ, a Q; = 1 ae ie Ip excess (84) 
NO 0 
TEN e 
AQwy oe 1 — anor I; excess ¢ (85) 
NO 


Aan »TBN 


ee erate 
vert eet 








Mic ys ™ 
YZ ZO 


(a) 





Fig. 22— The transistor in saturation. (a) Actual distribution of excess minority 
carriers. (b) Lumped approximation. The 7 and the T models use gy and gq; with 
lifetimes 7gy and 7g;; the Beaufoy-Sparkes model uses Qy sar and qgg = Agy + qr 
with lifetimes t—gy and rg, where rg is as defined in (88). 
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From this it follows that 


Qr 7 QAnoTcr 7 QIOT EI 7 AnoWaNn a So 
Since azo < 1, the charge-up ratio is somewhat larger than the ratio 
of the diffusion times of the normal and inverse transistor. 
The rate at which the two stores charge and discharge in saturation 
because of external step disturbances is described by the eigenfunction 
of the system 


eto Lt be | 1+ be), Lt bt Ba 


TBN TBI 


= 0. (87) 
TBNTBI 

If Byotz1/tezv >> 1, the two poles are far apart in frequency. Further- 
more, the high frequency pole contributes in most nonoscillatory cases 
little to the overall response. The higher pole or, alternatively, the s° 
term can then be neglected and a single time constant results described 
by 


_ G+ Bwo)rar + UL + Bio)ten _ Ten + Ter | 
ae Byo “1 Bro 1 — ayoero 


Using w-notation, one obtains the expression given by Ebers and Moll 


TS 


(88a) 


Ty = _ Wan FT War (88b) 


WeNnWe (1 rad AnoXr0) 


For large By and small 8; , 7s is approximately equal to 


Ty rail I a tax), (88c) 
TCI 

If tay K Tez, 1.e., if the carriers diffuse more easily from the emitter 

to the collector than vice versa, then the recombination rate 7,; on 

the collector side is mainly responsible for the overall decay of the 

excess base charge. 


6.5.1 Storage Time Calculations 


For first-order storage time calculations with the transistor driven 
into a steady-state saturation condition by means of an excess base 
current Jp excess , One can simplify the charge-control model to the one 
shown in Fig. 23. Storage time is the time it takes to deplete the store 
which is charged to a value of 


Qzes = Qr + AQy = I, excess? S — Iz excess Tay 1 OyoTer (89) 


l= Ano0@r0 
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Fig. 23 —Single-pole equivalent circuit for saturated transistor after Beaufoy- 
Sparkes. 


while, at the same time, this charge is exposed to an effective recombina- 
tion time of 7s as given in (88). In this form the model is identical 
with the classic Beaufoy-Sparkes model for the saturated transistor. 

In the general case one must refer to the complete model. 


6.6 The Beaufoy-Sparkes Charge-Control Model 


In the classic approach to charge-control theory, the starting point 
is, like in the diode case, the integration of the continuity equation (8). 
In comparison with the integration performed for the diode in Section 
2.3.1, the uppcr limit of integration has to be changed to x = W. The ex- 
pression 


i,(0, ) — 407, 9 = 22O 4 WO 
dt TBN 
obtained from the integration becomes that for the base current under 
normal, nonsaturated operation: 
_ dgqn(t) qv(t) , | 
tz(t) =i dt — Hea ’ (90) 
qv(t) is the total charge in the base. The next step being made is again 
the approximative assumption that the carriers redistribute themselves 
so quickly, that we can always assume steady-state distribution. 
(See also Section 2.4.) Mathematically, this means that in normal 
transistor operation both p(0,t) and z¢(¢) are proportional to the base 
charge qy(t). It can be seen from Fig. 19 that the same assumption 
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is Implied in the two charge-control models presented there, despite 
the fact that they were derived through entirely different procedures. 
(Instantaneous redistribution is, however, not implied in models that 
use more than one z or 7 section for representing the base.) 

The time constants are defined in the classic charge-control theory 
on the basis of the above-mentioned assumption of instantaneous 
carrier redistribution, i.e., in steady state . 


_ Qv _ Qv _ Qv 
Tan = Ton = oe TEN >= = (91a) 
and dynamically 
e d e » e 2 
tan >= dw + dw ’ ton >= dw ’ len = tay + ten . (91b) 
TBN dt TON 


The remaining three time constants can be defined likewise for the 
inverse transistor. Narud, et al’ have used such definitions in an equiv- 
alent circuit for the charge-to-current relations in the transistor. Beaufoy 
and Sparkes discussed this possibility in their original paper’ but chose 
to present two separate charge-control models, one for the normal 
active operation and one for saturation. In normal operation, the 
charge qy called ‘‘g,’’ is bounded by the value reached at the edge 
of saturation: 


ioe 


YB = On sat 3 where On sat — TBN « 
Byo 


In their saturated model, all excess charge which exceeds Qy sar 18 
lumped into one store rather than two; this charge ‘qs’ has a lifetime 
Ts = Qps/1p excess Which is identical with 75 as defined in (88). 

By lumping Aqy = qv — Qy sar and gq; into ggg , the Beaufoy and 
Sparkes arrangement provides only a minor short cut for calculating 
storage time, while sacrificing not only some of the physical under- 
standing, but also the possibility of mutual conversion with the other 
models. No relations have been given that would express the junction 
voltages in terms of the charges in the stores, and recourse must be 
taken to the Boltzmann equation to find expressions for the voltages. 

Throughout the literature the charge-control concept has been used 
primarily as a mathematical-physical tool. Extensions to higher-order 
effects are usually made by improving the simple continuity and 
transport equations stated in (8) and (9) and then carrying out the 
corresponding integration for the specific application. 
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VII. SOME REMARKS ABOUT THE EQUIVALENT CIRCUIT TYPE CHARGE— 
CONTROL APPROACH 


7.1 Use of the Charge-Control Models 


It is believed that the first-order approximation to a charge-control 
model in the form presented for the transistor in Fig. 19, combines 
the main advantages of the three basic approaches to modeling. The 7 
and the T-models are as easy to handle from an equivalent circuit 
point of view as the Ebers-Moll model. Instead of frequency dependent 
e’s and §’s, one additional current branch exists for each side of the 
transistor. Circuit problems are solved in the usual way by means 
of loop and node equations. The charges qy and gq; appear as circuit 
parameters which can either be calculated, if so desired, or else, elim- 
inated in the algebraic process. The store elements in the circuit are 
clearly defined by the circuit properties given in (22). 

The model provides all the features that have made the charge- 
control concept attractive in the past: quick estimates of switching 
times by Integrating the base current and equating with the charges 
needed to fill and deplete the stores. The general base current equations 
of charge control are directly read from Tig. 19(a) as 


Cx dv, 
dt 


Cre dv, 


T (92) 


: dqn d 
ee ae 


TBN dt TBI dt 


-{- 


Of course, there is no restriction to step inputs. The chore of cal- 
culating responses to a nonstep input is transformed through the model 
into a circuit problem. In complex cases the help of a computer will 
be required. 

Due to its direct relationship to the Linvill model, the charge-control 
model lends itself quite readily to extensions based on the physics 
of the device. This will be discussed in Section 7.3. 


7.2 Piecewise Linear Approximation of the Logarithmic Voltage Function 
The logarithmic voltage functions for the Junctions are of the form 
1 q/T | 
=e E =e 93 
x ON ” ve I,/C = yo&r0) ( ) 


For most practical cases, this can be approximated by piecewise linear 
functions, like in the diode case of Fig. 10. Except for small values 
of q/r, i.e., g/t not >>Io/(1 — e@yoaro) , one obtains 





(94) 
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Thus, the slope can be represented by a resistor r/\g, which may, 
like in Section 2.4.4, be taken as the average value 


dv OE et ere T 1 
a(2) ve Ag “7 BNO sie ere 
T 


where Jinax 18 the maximum forward junction current. 

It can be shown that in models which use the exponential relation- 
ship, the expressions of the form In(1 + 2) can be replaced by just 
In x if one simulates the majority carrier currents by special current 
sources as follows: 











(95) 


ab ins Ieo (leo) from internal base to collector 
1 — ayoayo 

and 
_1 = ao Teo (0) from internal base to emitter. 
1 — ayoezo 


This transformation is rigorous only at dc. However, in a piecewise 
linear analysis, as discussed above, the addition of one current source, 
namely J¢), becomes mandatory if the model is to be valid at very 
small collector currents. 


7.3 Hxtenstons of the Model 


7.3.1 Path Impedances, Leakage Resistors 


Like in the Linvill model, junctions and base material are clearly 
separated in the charge-control model. Therefore, it is a straight- 
forward procedure to add series path resistors, series inductances, or 
leakage resistances to models like the ones in Figs. 9 or 19. 


7.3.2 Junction Capacitors 


It has been indicated by the dashed lines in Figs. 9 and 19 how the 
junction capacitances are to be incorporated into the model. They are 
properties of the junction, but their currents flow as majority carrier 
currents through the bulk material. Hence, in Fig. 9, for example, they 
must be connected across the whole n-region. (Connecting directly 
across the voltage source would have no effect on the external prop- 
erties.) In Fig. 19 they lead to the internal base point. 


7.3.3 Higher Order than r-Transistor Models 


Another desirable expansion may be to replace the 7 structure of 
Fig. 19(a) by a double z or by some other higher-order approximation 
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to the original differential ‘‘transmission line’. This is of special im- 
portance, if emphasis is to be placed on charge redistributions in the 
base. By extending the model in such a way, the restricting assumption 
of instantaneous carrier redistribution is no longer implied. A qualitative 
example of an elaborate planar or mesa transistor model is given in 
Fig. 24, 


Cec 


i: DRIFT =F (Ve,Q;+Q2) LORRI St (Vesd2+4s) 


O01) 





Cpe OB Cac 


Fig. 24—Example of an elaborate high-frequency planar or mesa transistor 
equivalent circuit. 


7.3.4 Drift Fields 


If the charge-control model under consideration is being developed 
on the basis of physical phenomena such as in the model of Fig. 19(a), 
the contribution from drift effects may be represented in the same way 
as has been proposed by Linvill [Ref. 7, Sections 2, 3]. As a direct con- 
sequence of the transport equation (2), drift can be represented by a 
current source added in parallel to the diffusion current source. In terms 
of the r-g-c transmission line representation, discussed earlier, drift con- 
sideration amounts to a resistor in parallel with the series diffusance 
resistor r. This was used in a recent paper by Bloodworth.” 

Alternative methods of representing drift effects in conjunction with 
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conductivity modulation are presently being investigated; results will 
be published later. 


7.3.5 Base Width Modulation 


Base width modulation can be taken into account by replacing the 
basewidth, especially the collector section W., by an expression 
W.(1 + A), where A is some function of the junction voltage. Equa- 
tions (52) and (53) show the dependence of the branch currents on W, , 
from which we can readily derive the required modification of the 
charge-control model in Fig. 19(a). 


7.3.6 Multtple-Layer Devices, Multiple Storage 


In accordance with Linvill’s proposal, storage in more than one 
region can be simulated by considering that the minority carrier cur- 
rent on one side of a junction becomes the majority carrier current 
in the adjacent region. An example is shown in Fig. 25. This figure 
represents the charge-control model for an npnp device. Avalanche 





L ELECTRON 
= — ras \ 
Vy,=f (pn) L=f (qpn,4pr,M) Vo=f(qp1) \ OGATE 


© | 
/ 


EQUAL 


+) Ct) am \ 
/ 


/ 


/ 
/ 
/ 
/ 
: (v) > (vA 
lL HOLE / 





‘GATE Vyo=F (Qn1) L=f (qnn,4nt™) Vyg=F (Qnn) 
ee ery 4 
JUNCTION J; D-REGION JUNCTION Jo N-REGION JUNCTION J3 


Fig. 25—— Charge-control model for npnp device. The model for an npn-transistor 
with storage in the collector can be obtained from this model by omitting the part 
to the right of the dashed or the dotted line. 
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multiplication may be considered by adding a multiplication factor M7 
to all hole and electron currents flowing through the junction of interest 
(usually the center junction J,), as indicated in Fig. 25. M7 is a function 
of the voltage across the junction. 

By omitting the last electrode, an npn transistor with charge storage 
in the collector is obtained. 


7.4 Iistablishment of a Large-Signal Model 


The question naturally arises as to how one arrives at a numerical 
model. There is no clear-cut answer to this question, since the procedure 
to be taken depends on whether the informations available are pre- 
dominantly physical or electrical in nature, whether a computer is 
available or not, etc. The following outline can, therefore, only be 
considered as a typical example. 


(z) Obtain de measurements which yield information on junction 
characteristics and electrode resistances. All measurements must be 
made under widely differing drive and load conditions. 

(ic) Add information from device manufacturer to establish first- 
order de model. (If necessary, convert to Linvill model.) 

(221) Add dynamic parameters, such as capacitances, as far as they 
are known and establish first-order dynamic model. 

(iv) Use computer to improve numerical parameter values by 
matching frequency response curves or switching data in the active 
region with the model. 

(v) Use computer to match large-signal nonlinear switching data. 

(vt) Check model with switching measurements under different con- 
ditions, such as extremely low, extremely high and medium input and 
output impedance levels for various drive conditions. Improve model 
basically and numerically as necessary. 

For purposes of device design, more emphasis is generally placed on 
the simulation of higher-order effects than in model building for circuit 
design where, especially in the case of integrated circuits, it 1s necessary 
to trade accuracy for simplicity. 


VIII. CONCLUSIONS 


The differential Linvill model stands out among all models as the 
most perfect one. Whereas the lumped Linvill model is the most suitable 
model for the device physicist, the circuit engineer usually prefers a 
more circuit oriented approach. It is felt that the charge-control equiv- 
alent circuit approach is well suited to combine the main advantages 
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of the various models: It is as easy to handle as the Ebers-Moll model, 
yet bears the close relationship to the physical phenomena of the 
device inherent in the Linvill model. It can also be extended easily 
to include higher-order physical phenomena. 

Despite the difference in basic philosophy underlying the creation 
of each of the three classic modeling concepts (such as lumping, super- 
position, and integration), they are equivalent with respect to their 
current relations and to all de properties. When compared at the same 
level of complexity, equivalency with respect to time dependency of 
the junction voltages exists between the two charge-control models 
and the Linvill model, but not between these models and the Ebers- 
Moll model. | 

In the Ebers-Moll model the effect which storage exercises on voltage 
cannot be included. The hydrib use of both time and frequency domains 
in the model may also be felt as a disadvantage in some applications. 

At the first-order level of approximation, the charge-control equiv- 
alent circuit can be converted into the Ebers-Moll model, the Beaufoy- 
Sparkes model, and into the Linvill model (in the latter case the base 
volume is a constant which must also be known). Thus, the charge- 
control model serves as a bridge between the various models. This can 
be very useful in establishing a model, since both physical and electrical 
information can be incorporated easily into the model. 

The diode charge-control model has been found very useful for 
analyzing storage effects in diodes. 

Because of the close relationship to the physical phenomena in the 
device, extensions to larger complexity can readily be accomplished. 
We may interpret the charge control equivalent circuit as simply a 
circuit-oriented form of the Linvill model. The basic ideas and pro- 
cedures that are used in converting diode and transistor linear models 
into equivalent charge-control models can be applied to many other 
semiconductor devices. 


APPENDIX A 

Switching Time Calculations for Ideal Charge-Storage-Step-Recovery 
Diodes 

(Example for use of charge-control model) 

A.1 Equivalent Circuit (See Fig. 26) 

A.2 Generator Source Current (See Fig. 27) 

A.3 Diode Model (See Fig. 28) 
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L(t) V_(t) 


Lo(t) + 
Lo(t) Ro Rr RLSfv 


— _ 


Fig. 26-— Equivalent circuit for charge-control model. 


Lo(t) 





Fig. 27 — Generator source current. 


eee L(t) 


q/7 
SWITCH CLOSES AT t=O 
OPENS WHEN q=0 


Fig. 28 — Diode model. 


A.4 Forward Operation 


1) = 22 + sq@ = =. 


si 
From this follows 
I; 


+7) 
qt) = rIy[1 — exp (—1/7)] 


Q(t,) ad rf ,(l — €xp (—1,/7)] 
A.5 Reverse Operation 


Qs) 


For simplicity of writing, t = ¢, will now be referred to ast = 0: 
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1) = 29 4 sag — a) = -®. 


T S 


From this follows 


Q(s) a sQ(t,) — Ip 


fe 


qt) = Q(t) exp (—t/7) — rLp[l — exp (—¢/7)]. 


Step recovery occurs att = T,, when g = 0 


exp (—T,/7)[Q(Z,) + tle] = rlre 
T= 7 In }1 + 22 | - 7 In G+2n — exp (— /n). (96) 


A.6 Graphical Representation (See Fig. 29) 


q(t) 





Fig. 29— Graphical representation. 


A.7 Special Cases 





(1) tp: QC) = Ipr 
Pe ah (1 4 Z| (97) 
WW. ES Be ahs (98) 


(ut) Ip>>Tp and th 0: T, ==>. (99) 
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APPENDIX B 


Comparison of Notations 


Table I lists comparisons of the notations used in this article. The 
first column lists the notations used in this article while column A 
lists those used by Beaufoy and Sparkes.” Column B lists the notations 
used in Physical Electronics and Circuit Models by P. E. Gray, et al;”* 
SEEC Series, 2. 


TABLE I—CoMPARISON OF NOTATIONS 


A B 
This paper Beaufoy-Sparkes SEEC 
Qn ge(for gy < Qnsar only) QF 
qr dr 
Agn + qr 
QBS 
where Aqy = Gy — Qysar 
TBN T's TRF 
TON Tc TF 
1 1 
TEN Tr 1 / (2. ae 4) 
TBF TF 
TBI TBR 
1 1 
a Acme 
TBR TR 
TEI TR 
TS T's TSL 


LIST OF SYMBOLS 


Lower-case letters are used for time variables, capital letters are used 
for steady-state values or Laplace transforms of values. 


A, Ai, Ae cross-sectional areas 

Ai, Ay, Aor, Ace four-pole parameters 

a, b, c, K constants 

Cf analog capacitance; same per unit length 

Cre, Cre emitter and collector junction capacitance, 
respectively 

D, hole diffusion constant 

D, electron diffusion constant 

e magnitude of electronic charge 

E electric field intensity 


analog shunt conductance; same per unit 
length 


—) 
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Beast 
H. 

Hi 

1/H. 


41, 42,43, U4, U5 





Up 

Lp excess ) I; excess 
te, Le 

tg, Le 


Ip, Lr 


by 

ty 

L's 

U6 ? I¢é ’ Un Ie 

lor ? Lor ) lar 5) Lup 

Lo, Leo, Leo, Lero, Lnro 
Ie Sat 


li. 5) Lita 5) Us 
tpn » tEN } lcN 


In 
Ip 


L, 


Mm, M4, , Mz 


n 

Np 

Np» Npo 

p, P(s) 
Pn 

Po» Pno 
q,Q 

G1» W255 Ym 
WN 5 Qy 

Qi ; Qr 


Or 
~] 
ow) 


lumped combinances 
combinance per unit length 
lumped diffusance 
reciprocal of diffusance, per unit length 
network branch currents 
base current 
excess base current in saturation 
collector current 
emitter current 
forward and reverse diode switching current, 
respectively 
electron current 
hole current 
diode saturation current 
ae branch currents as defined in 
Figs. 17(a) and 17(b) 
de junction saturation currents 
collector current in saturation 
currents through combinance, diffusance 
and storance, respectively 
base, emitter, and collector current in normal 
transistor operation 
electron current density 
hole current density 
Boltzmann constant 
diffusion length for holes 
constants relating lumped carrier density 
to carrier density at junction boundary 
electron density; excess electron density 
excess electron density in p-region 
values of n and n, in thermal equilibrium 
hole density or excess hole density 
excess hole density in n-region 
value of p and p, in thermal equilibrium 
charge 
lumped charges in base region 
charge in normal store 
charge in inverse store 


Aqv , AQy, Agr, AQ, Qzs additional charges stored due to saturation 


Qy SAT 


limiting value of Qy, reached at edge of 
saturation 
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Gio» Ymo 
r;7r 


QN ; Qy 


Ono; “ro 


By , Br 


Byo ) Bro 


Y 
r 


Bn, Up 


, 
p> Tay Tbh» T12 


Ts 

7, = TBn 
To = Ter 
TON » TEI 
TEN » TCI 

WaNn ) War 
Wen » Wl 


total minority carrier charge in equilibrium 
analog resistor, same per unit length 

small signal junction resistance 

Laplace operator 

stores = storances 

storance per unit length 

time 

absolute temperature 

voltage 

externally applied junction voltage (ex- 
cluding resistive drops) 

collector and emitter junction voltages 
lengths denoting sections in neutral region 
base width 

neutral region length variable 
characteristic impedance 

normal and inverse ac current gain in 
common-base connection 

de values of ay and a; 

normal and inverse ac current gain in 
common-emitter connection 

de values of By and £, 

transmission line propagation constant 
short for e/kT 

electron and hole mobility, respectively 
recombination time constant in p-region 
diffusion time constants 

approximative effective recombination time 
constant for excess charge in saturation 
recombination time in base under normal 
operation 

recombination time in base under inverse 
operation 

normal collector and inverse emitter diffu- 
sion time constants, respectively 

normal emitter time constant (= 1/w,y,) and 
inverse collector time constant (= 1/w.z), 
respectively 

common-base angular cut-off frequencies 
common-emitter angular cut-off frequencies 
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Generalized Optimum Receivers of 
Gaussian Signals 


By T. T. KADOTA 
(Manuscript received October 28, 1966) 


Optimum reception of two zero-mean Gaussian stgnals 1s accomplished 
by comparing a quadratic form [f x(s)H(s,t)x(t) ds dt in the observable 
waveform x(t) with a predetermined threshold, if the symmetric kernel 
H(s,t) can be given as a square-integrable solution of 


| Ris wWH(u,v)R.(0,) du dv = R,(s,t) — R,6,)), 


where F,(s,t) and R.(s,t) are the covariances of the two signals. In this 
paper, we generalize this result so that oi. ff x°?(s)Hin(s,t)x”” (t) ds dt 
as the quadratic form to be used and {H,,,(s,t)} ts given as a formal solu- 
tion of 


L m 
X [f Se Reptntus) & Rab) du dv = Ralo,t) — RyG). 
L.m 


In other words, the generalized quadratic form 1s in the derivatives of x(t) 
as well as x(t) ttself and the kernels H,,,(s,t) constst of two-dimensional 
6-funcitons in addition to square-integrable functions. This result 1s ex- 
tended to the case of two nonzero-mean signals and then to the case of M 
Gaussian signals in noise. 


I. INTRODUCTION 


Consider the problem of discriminating between two zero-mean 
Gaussian signals by observing the sample function z(¢), 0 S$ ¢ S 1. 
We assume that their covariances A,(s,t) and F,(s,t) are continuous 
and positive-definite on [0,1] X [0,1]. According to previous results,” 
if the integral equation 


ih [ Ri(s,wH(uv)R20,1) du dv = R,(s,t) — Ris), OSstS1, (1) 
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has a symmetric and square-integrable solution H(s,t), then the fol- 
lowing decision scheme minimizes the error probability: 


1 1 
choose f,(s,/)_ if | | x(s)H(s,t)x(t) ds dt < ¢, (2) 
0 0 
choose f,(s,t) otherwise, 
where 
c= 2log— — >, log d; , (3) 
Qs +=0 


in which a, and a, are the a priort probabilities associated with the 
two signals, and A; > 0,7 = 0, 1, 2, --- , are the eigenvalues of an 
operator R7?R,R>?.* 

Unfortunately, existence of a square-integrable solution of (1) is 
too restrictive a condition. Thus, relaxation of the condition, which 
amounts to generalization of the quadratic form of (2), is desirable. 
In this paper, we accomplish this in two ways: one is to allow H(s,t) 
to contain 6-functions as well as square-integrable functions, resulting 
in the generalization of the structure of the quadratic form; the other 
is to consider the derivatives of x(t) as well as x(¢) itself, thus generalizing 
the elements of the quadratic form. The result is extended to the case 
where the means of the two signals are nonzero, and is further extended 
to the case of 7 Gaussian signals 1n noise. 


II. GENERALIZED OPTIMUM RECEIVER OF TWO ZERO-MEAN GAUSSIAN 
SIGNALS 


Consider the following generalization of the quadratic form of (2): 
r 1 1 
a@ = DO ff 2? @mals,a”( ds at, (4) 
Z,m=0 40 0 
where x‘ (2) is the [th derivative of x(#), and 


H1,AS,t) = d Qjkim d(s == S;) b(é = Sz) 
1,k=1 


* More precisely, A:, @ = 0,1,2,..., are the eigenvalues of the extension 
of R|?R.R;,* to the whole of £2, where R, and FR, denote the integral operators with 
the kernels F,(s,t) and R2(s,t), and £2 the space of all square-integrable functions on 
[0,1]. We recall that existence of a symmetric, square-integrable solution of (1) 
implies that Ry?R.RT* has a unique bounded extension to the whole of £, having 
eigenvalues {\;} such that 0 < IlP_p Aa << o&. 
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. pS Cee ee Oe ROR Cee 
+ him\s) 5(s _ t) + H,,,(s,2), (5) 


in which a;,;,, are real constants, and OSs, , ¢; <1, and hj,,,(t), hj1,,(t) and 
him(t) are square-integrable functions on [0, 1] while #,,,(s, £) are square- 
integrable functions on [0, 1] < [0, 1]. In writing (4), we have assumed 
that almost all sample functions of both signals have rth derivatives.* 
Note that the nonsquare-integrable part of H,,,(s,f) consists of three 
types of two-dimensional 6-functions: (2) those at points and their 
mirror images with respect to the diagonal s = J, (2) those along 
horizontal lines (¢ = constant) and their mirror images (s = constant), 
and (777) those along the diagonal. By formally substituting (5) into (4), 
we obtain an explicit form of Q(z), namely, 


Q(x) = > LS Ajkimet : ‘(s; a” (S,) 


t,m=0 


+ F 2) [tran + FrrnlD](D at 


+ fo Ohne at + ff x? OA rmlse( ds at} (6) 


As the corresponding generalization of the integral equation (1), 
we consider the following: 


~ f [2 Ai oes (,U) A in(u Oe = Ral t) du dv 


L,m=0 
= R,(s,t) — R,(s,d), Ossil. (7) 
Again, through formal substitution of (5), (7) becomes 


r n1 a! a” Ne 1 
>» >) Qikim at! POS D) lice au” R.(s,) S=8k + 2, i 


1,m=0 





[i ih i(S, t) |.— ~T, hsam(t) Ratu, t) 


+< aut +R i(s whijrnltl) a Reals, t)|.~ | au 


* A simple sufficient condition for this is existence of (d°°t?/ds™tat"+!) R,(s,t), 
4 = 1, 2.4 
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+ f Zr oaheny racy aut ff 2 Rea fratus 


< Ro(oyt) du iv} = F,(s,t) — Ri(s,0), 0Osst<1, (8) 
where we have assumed that (0°’/0s’0t")R,(s,t) and (0°"/ds"dt’)R,.(s,t) 
exist and are continuous on [0,1] X [0,1]. 

Unlike H(s,t) of (2), which is uniquely given as the symmetric, 
square-integrable solution of (1),* the defining elements of Q(x) (..e., 
{Gizim}, {Si}, {ti}, (hiom(}, (Rirm()}, {hin(t)}, {Him(s,é)}) cannot be 
uniquely determined by (8) in general for a given pair of covariances 
R,(s,t) and &.(s,t). Nevertheless, we can establish the following: 


If (¢) R,(s,t) and R,(s,t) are positive-definite, 

(iz) (0°"/ds' dt") R,(s,t) and (0°"/ds"dt’)R.(s,t) are continuous, 

(222) for almost all sample functions both signals have rth derivatives, f 
and 

(iv) there exist some set of finite sequences {a;x1n}, {S;}, {t;}, {Rio}, 
{Airm(t)}, Arm(t)} and {H,,,(s,t)} which satisfy (8), then the decision 
scheme (2) with fofo x(s)H(s,t)a(t)dsdt replaced by Q(x) of (6) is 
optimum. 


The proof is based on two measure theoretical facts: (z) two prob- 
ability measures P, and P, corresponding to two Gaussian signals are 
either equivalent or singular,{°’? and (72) if they are equivalent then 
there is a special random variable called the Radon-Nikodym derivative 
(dP./dP,)(x), in terms of which the optimum decision scheme is 
specified as follows:" 





on dP, Oy 
choose £&,(s,f) if aP, (a =< . 


choose £&#,(s,t) otherwise. 


Hence, in the Appendix, we first prove that existence of {Ajzim}, (8;}, 
{ti}, (hirm(Q)}, (Ajim(t)}, {Ain(é)} and {#;,,(s,t) satisfying (8) implies 
equivalence of P, and P, . Then, it follows that the eigenvalues \; , 7 = 





* The uniqueness of H(s,t) follows from positive-definiteness of R;(s,t), 7 = 1,2, 
and square integrability of H(s,t). 

+ Continuity of (6?"/ds"at") R;(s,t), 7 = 1, 2, and existence of 2‘(¢) for almost all 
z(t) may be replaced by a simpler but stronger condition that (02"*?/ds™*1at7t1) 
R,(s,t), ¢ = 1, 2, exist. 

{ From the communication theoretical point of view, singularity corresponds to 
the case of “perfect reception”? where error probability vanishes. For the mathe- 
matical definition, see Ref. 7. : 
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0,1, 2, --- , exist.°’* Next, we explicitly obtain \,; from (8) and show that 
0<[[2., \;<. Thus, the threshold c of (3) is well defined. Lastly, 
_we prove that 





oe) = (Tx) ew 40 © 


for almost all z(t) of both signals. Then, by substituting (9) into the 
above decision scheme and taking the logarithm of both sides, the 
assertion is immediately proved. 


IiI. EXTENSION TO TWO NONZERO—-MEAN GAUSSIAN SIGNALS 


The preceding result can be extended to the case where the means 
of the two Gaussian signals are no longer zero.* Let P,, and P22 be 
two probability measures corresponding to two Gaussian signals with 
means ™m,(t), m.(t), 0 S ¢ S 1, and covariances R,(s,t), R2(s,t). m,(2) 
and m.(t) are assumed square-integrable while the assumptions on 
FR,(s,t) and F.(s,t) remain the same. Introduce a third measure P., 
corresponding to a Gaussian signal with mean m.(¢) and covariance 
F,(s,t). Then, P,, and P2. are equivalent and 


ope) = Sp) 2 © 





for almost all x(¢) of all three signals, if and only if Ps. is equivalent 
to P.,, which in turn is equivalent to P,,. According to a previous 
result,’ if there exist finite sequences of real numbers {d;,}, and {é;}, 
0 < i; S 1, and square-integrable functions {g,} which satisfy 


r na a: 1 go! . 
De b> Git Bet Basi) leas = | — Ri(s,t)Gi(s) as 
j=1 8 0 Os 


1=0 


for almost all x(t) of the two signals, then P,, and P.; are equivalent 
and (dP.,/dP,,)(«) = exp [L(2)] for almost all x(t) of the two signals, 
where | 


L(2) = 01S an E () - mde) tml) | 


| [ o at) — mall) + mall) Io, cp it. (11) 7 


* This extension follows the development in Ref. 3, pp. 1628-1629 and pp. 1636-— 
1637. 


t This is the “sure signals-in-noise”’ counterpart of the result mn Section II, 
namely, the generalized optimum receiver of two sure signals in Gaussian noise. 
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The remaining half of showing the equivalence of P., and P.. and 
obtaining (dP2./dP.2,)(x) is accomplished simply by replacing x(t) with 
x(t) — m.(t) in the result in Section II. Thus, upon combination, we 
conclude that, if there exist a set of finite sequences {4;,}, {t;}, {G:(t)} 
satisfying (10) and another set of sequences {d;x1m}, {8}, {t;}, {Rim@}, 
{Riim(t)}, {Rim(t)}, {Hin(s,t)} satisfying (8), then the optimum decision 
scheme for this case is specified as follows: 


choose m,(t), Ri(s,t) if 2L(c) + Q(z — me) <.¢, 


choose mz,(é), R2(s,t) otherwise. 


IV. EXTENSION TO M GAUSSIAN SIGNALS IN NOISE 


The above result can be further extended to the problem of dis- 
criminating among Jf Gaussian signals in Gaussian noise.* Let m;(¢), 
R;(s,t) and a;,7 = 1,2, --- , M, be the means, covariances and a priori 
probabilities of the signals, and R,(s,t) the noise covariance where 
the noise mean is assumed zero. The assumptions concerning m,(t), 
R,(s,t) and #,(s,t) are the same as in Section III.t Denote by P;; the 
probability measure corresponding to the 7th signal plus the noise, 
and by P, the measure corresponding to the noise alone. Then, according 
to the theory of the generalized maximum likelihood test," if each 
P,;; is equivalent to P, ,t then the optimum decision is to choose that 
m,(¢) and R,(s,t) for which a,(dP;;/dP,)(x) is maximum as a function 
of 7.8 Observe that, if the 7th signal plus the noise and the noise alone 
are interpreted as the two Gaussian signals of Section III with means 
m,(t) and zero, and covariances F,(s,t) -+- &;(s,t) and #,(s,t), then 
the condition for equivalence of P;; and P, and the expression for 
(dP;;/dP,)(x) are obtained simply by the following changes: m,(t) — 0, 
m2(t) — m,(t), Ri(s,t) — R,(s,t), R.(s,t) — R,(s,t) + R;(s,t). Thus, 
we conclude that, if for each z there exist a set of finite sequences 
{Gar}, {é ti} and {Gir(t)} satisfying 


r Nis F Q! 1 a’ 7 
2 » Qijl jaf olSst) lames = [ pet Peolsst) Gurls) is| = m,(t), 


t=0 


0 


IIA 
IIA 
La 


* This extension follows the development in Ref. 10, pp. 2192-2194. 

1 &;(s,¢) need not be strictly positive-definite. 

sf Equivalence of P;; and P, corresponds to the condition that the 7th Gaussian 
signal cannot be detected perfectly in the presence of this noise. 

§ If a(dP;;/dP.) (x) becomes maximum at more than one value of 7, choose the 
lowest of such z-values. See Ref. 11. 
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and another set of finite _ Sequences {Asszim}, {Si}, {ti}, {hein}, 
{Aisim(6)}, {Air (t)} and {H;,,,(s,t)} satisfying 


r nia Q! a” 
ay 7 Qijkim apt BeolSst) ewes ij as™ [Fe (s, t) = R; i(s, b) eaves ik 


i,m= 


+¥f E Roberson) = (Ry(u,t) + Bile.) 

a 2 ; Rls hiss nll) - he (s,t) + B,(8,t)) |e: Jax 
+f[2 2p, (s,ui) hea (t) 

co [R.(s,t) + R.(s,t)] du + | | I | 2, Rs.) Harn(ua) 


ay [R.(v,t) + R,(v,d]| du awh = F,(s,b), 0<s,¢ = 1, 
b 


then the optimum decision is to choose that signal (m;(t), R;(s,t)) for 
which 2L,(x) + Q;(@ — m;) + ¢; is maximum as a function of 2, where 
L(x) and Q;(x) are defined by (11) and (6) with dj,, é;, 9,(é), and 
Qikimy S33 o hiim(8), hiim(t) him(2), H,,.(8, t) replaced by a Qij1, t sis Gilt), and 
Qijkim) Siz, b a7) hijim(), hiitm(b hirm(t); A iin(s,t), respectively, and 


c¢; = 2 loga; — > log An”, 
n=0 
where \{”, n = 0, 1, 2, --- , are the eigenvalues of the extension of 
I + R;?R;R;?* to the whole of £&, . 


APPENDIX 


Let P, and P,. be two Gaussian measures associated with a separable 
and measurable process {x(t), 0 S t S 1} with means zero and co- 
variances F,(s,t) and &.(s,¢). 


Theorem: Suppose R,(s,t) and R.(s,t) are (strictly) positive-definite, 
(0°"/ds’0t")R,(s,t) and (d°"/ds"dt’)R2(s,t), O S r < ©, exist and are 
continuous on [0,1] X [0,1], and almost all sample functions have the 
rth derwaties with respect to P, and P,. If there exast a set of finite se- 
quences* {Ajxim} {8;}, {ti}, (hirm(O)}, (hirm(} {hin(} and {H1n(s,t)} 
which satisfy (8), then 


* The definitions of these sequences are given after (5). 
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(2) P, and P, are equivalent, 2.e., P; = P2, 
(iz) (10) holds a.s., [P,, P2].* 


Proof: For simplicity, we introduce the following notations: 


lg a” 
R;.:(u,t) = pat Pil8st) emu ’ Rs 4m(s,v) a7, ajm Fea(8,t) lens ’ 


L+m 


CO ae rye Tae RCD cniy 21g? 


Tr ni 


K,(s,b) an > > Askimderet(8,8;) Ro_n(S; ’ b), 


i,m=0 7,k=1 


K,(s,l) = 3 3 [Riet(s,t;)(Reemhsam) (E) + (Ricthiim)()Rom(t; , I], 


t,m=0 


K,(s,t) = = [ RewncoR.ca pay. 


r 1 1 

Rens / i R,,1(8,) LE n(tt 0) Raen o,f) dea do. 
l,m=0 “70 0 

Note K,(s,t), 7 = 1, 2, 3, 4, are square-integrable. Again, we delete 

the arguments s and ¢ of the kernels to denote the corresponding integral 

operators. Thus, (8) becomes 


4 
a K, = Bs = R, ’ (12) 
t=1 
hence, 
4 
RRR; — I = >> Ri'K,R7. (13) 
i=l 


(1) To establish P, = P, , it suffices to prove that R7?R.R{? is densely 
defined on £ and R7?R,R{? — I is of Hilbert-Schmidt type, i.e., 


||RUIRRy? — I|| < 08? ?F The principal tool to be used for this 
proof is the following expansion: 
Risin) = DawfPOnk’O, OSlmsr, (14) 


uniformly on [0,1] X [0,1], where nu; > O and f;(@),7 = 0,1,2---, 
are the eigenvalues and orthonormalized eigenfunctions of FR, . 
To prove that R7?R,R{? is densely defined on £, , it sufficies to show 


* “a's.[P1, P2]” is the abbreviation of “almost surely with respect to Pi and P2’’. 
+ Hdl ae. the Hilbert-Schmidt norm of an operator A. 
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that R7?R. is bounded since R7? is densely defined. Now, by applying 


the formula ||A||? = tr A*A = 50; (f;, A*Af,) to the individual terms 
of R7?K, first, we obtain 


. r ny 
[| Rr*K, || Ss De = | Gina | 
i,m=0 7,k=1 


ip | [ (Ri*f,) (u) By. (u,8;) 





2 1 1 
i Rogm(S, ; U) Roim(U,S) an | 
0 


r n1 


= DD > | Qixem | | Ds Bs | f(s) |” Roemen(S, » S,) |? 
l,m=0 j,k=1 t 
r ny : 
ona 2 >. | Qiklm | | Bigizi Ge8,) R35 .mealS; ) Sx) ie 
i,m=0 7,k=1 


where (14) is used for the last two equalities. Similarly, 


| RK, | < p> > { (Rist (t; ’ t;) (itm eee ree a 
»m=0 j=l 


+ [Nj mt Ae Per ee cere OF ’ pie 
| RI?K; | < p> | tr (Ravi x asics) P. 


\| Ri?K, || S y= Ge CR piste dee weal 
»m=0 


where Ry,.:,:,,(8,t) = Rim (S)Ryo1,1(8,t)Rin (2), Ristitm = HmtRytitHin - 
Hence, from (13), ||Ry7?R.|| << ©. 

To prove ||R{?R.R;? — I|| < ©, we apply the formula ||A||? = 
>= ||As.|/? to the individual terms of R7?K,R7? first. Thus, we obtain 


I Rr" || SDE DE | aieem | | D2 || Re Roem sedeifeP(6,) UIP 
»m= 7,k5 2 


Tt Ni 


= » > | Qikim | | Ry .t.0(s; , 8;) b | Hs lospal5ey) II, 


l,m=0 7,k=1 


where (14) is used twice, and R7?R2,m(-,s,) denotes the result of Ry? 
acting on an s-function R.,»(s,s,). By differentiating both sides of (8) 
with respect to ¢, we obtain 


Ro m(S,8,) re Ri im(S,Sz) + p> Qjeimler,t(S,8;) l,m ym (S, » Sx) 
»m=0 


»m=0 7=1 


a D> > [Fy 12(8,t;)(Rosmemlsim) (Sk) + (Ri ithjim)(8)Rosmim(t; » Sx) 
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r 1 
“i 2 il Ry its whim (Ut) Rom ,m(UyS:) du 
i,m=0 YO 


r 1 1 
+ do ff Richa) Henltet)Ravmen(oss) du do, (15) 
L,m=0 70 0 


Thus, 
| Ry*Ro,m(+ Sz) | S | Fer amsm(Sp , Sz) |? 


r ni 
“ir pz) a | Qikim | | Ris; , §;) |? Rosmen(S; , Sx) 
0 7,k=1 


+ p> ae | Raster(ts 5 #4) |? | (Reemenljrm) (Sx) | 


»m=0 7= 


-{- ii. Rrgtethiim)? | Roemen(S; , Sx) |} 
+ DT Beemem(e +) Bretetm Roamen(: 81)? 


ae 3 (Rosmim(Sz y *), Risteig Rogmim(: S,))* 


L,m=0 
< 0, (16) 
Hence, 
|| RRy? || < &. 
Similarly, 


|| Ri*K.R;* || S > 2 [| Raster(ts st) |* || Rr *Reemhjrm |I 
+ (eis sighed) | Ry Rarn(- yt;) [{]. 
From (15), 
| Ri Rejolqx | < (hitm Rigarahei)? 


as x | Qjeim | | Rissee(s; 5 8) |? | (Rosmembsim) (Sx) | 


l,m=0 


+ = pay Riidlhst) || Craw hiwehaay| 


li,m=0 j= 


= (iim ’ Ristahiis)? | (Ro,m:ml;1m)(t;) 1] 
st 2s (him ’ ane eee Rhian) 


a 2 (him Rew chien) 


t,m=0 


< 0, 
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and, from (16) 
|| Rr*Roe(- 3) || < ©. 
Hence, 
|| RetKeR,* ||-< eo, 
Similarly, 
|| RrKoR;? || 

















2 > br Rf Raa(- than fe (a) au | 
t,m=0 t 0 
SE {le Bare Rane) PE SS Lane 


| Ry gt’ yt’(s; ’ S8;)(Rogm'ym(S, ’ 5 Retin Raogmym'(+ Sy) |? 
+ > >> [| Ryu "eb? (CL; ’ E)\(Bogmym? hjv m’ 9 Ristet.m Rogmyem' Retr) |? 
»m’=0 j=1 


+ i. isha) eieulcel G ) ays Brace, awe) 


We > | tr Citic as Re setyal agiew Roanyn’) 3 


~m'=0 
; 
2 





g 
o~w a 
si > | tr Ler et etm? Rogm’ pmbtygtet, m Pagmym') 
Ll’ ,m’=0 


< 0, 
Similarly, 
|| ROKR || 


p>) > By | Ry Rem mifs 


IIA 





r 


{ tr Ristil cs Ri, 


A 
M 


a > oy | Qjaa’ m! | 
t,m=0 '=0 7, 


>| Rustrar(s; » S;)(Roqm'pm(S, ’ 2) Miatyl x Rogmym'(* Sz) |? 





le » Dy [| Rist'r(t; ’ bi) item’ ce actaltde petgtics Ragmym'hj 1m’) |? 
L';m'=0 fal 
- | Rigrrms Rast et Rjrrm) (Rosman tj ’ “) Rigiytcss Ry,mym'(> ,t;)) |] 


a ; 

ie A owe 1 

#15 » | tr (Lor 9t' et’, m! Rogm'pmbtystet,m Rogmim’) |? 
l’,m’=0 
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r 
ome —~ 1 
+- Ds | tr (Feit et sm Lis cmt do igig® ben Ro, nim’) | 
l’,m’=0 


< om, 
Therefore, from (138), 
(| Res — fF || <= x: 


(it) We have established in (2) that R7?R,R{? is bounded and densely 
defined on £,. Hence, it has a unique extension to the whole of £2, 
which we denote by M/. Since 4 — TI is a Hilbert-Schmidt operator, 
M has eigenvalues and orthonormal eigenfunctions, which we denote 
by A; and ¢;(é), 7 = 0, 1, 2, --- . Note 0 < A; S | M |, where | M | is 
the norm of 17. Then” 


Rails) = » (Rig) (9 (Rig). 


Osl,mSsr, (17) 
Rostn(s,t) = DY r(Rig,) (9 (Rieg) , 


uniformly on [0,1] * [0,1]. 

Let {¢,,} be sequences of functions in the domain of Ry? such that 
y; = lim. ¢;, for each 7. Multiply both sides of (12) by (Ry*¢;,)(s) 
and (Rj%¢,,)(t), integrate with respect to s and ¢, and let n > o. 
Then, the four terms on the left-hand side become 


(Ry pin ’ Ki Rye) 


r ny 


» >, tithe os »Rrier(: 8;)) (ogm(S;, ’ )s Ry ein) 


t,m=0 j,k=1 


l 


r n1 


= Dd Do Onin pe (vin 5 9») (Rig,) (s;) Zs r (Rig, (Sie » Pins 


i,m=0 7,k=1 


where (17) is used for the second equality. By virtue of (17) again, 
we can define an s-function R?,:(s,u) for any u e [0,1] by 


R?,(s,u) = lim. >> ¢,(s)(Rig,) (u). 


n> v=0 
Then 
zx (Cin ? v,) (Rig,) (s;) = (Sin ’ R? ,(- 8;)) 5 


a d,(Rig,) ™isd(¢, Pin) = (R} ,m (Sy ’ -), Mein), 
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and 
lim. (ys , Ri(+ 8) = @ , Ris) = (Rig) ,), 
mu (Ri.m(Se , +), Mein) = (Ris » +), Me) = (Rigs) (sd). 
Hence, 


rT 


lim (Ryton , KRU em) = A pa] >» Ajxim(Rig;)“ (s;)(Rig.)” (s:). 
m= k=l 


NCO 


Similarly, 
lim (RY*¢in , KoRi*gin) 


n—?0Q 


= lim > 3 [(R7? Pin | Rul: t;))(Ropmhjrm ’ Ry +0 in) 


nc l,m=0 j= 


+ (Ry vin Riaz) Rae: ’ ), Ry tin) 
= =; >» pa (Ri Pi OE, ) (Rig, ae iim + hirm) 


lim (Ry49;, , KsRi*ein) 


n-o4 


= Fim fib ein y Rist) ham(t aunty), RM) de 


n—0 L,m=0 


=D fled hina Rie) Co) du 


L 


lim (Ry Gin ’ K Ri gin) 


naw 


r 1 1 
=lim »> ime) (Ryteim 5 Rast(-,u))Him(uv)(Rosn(v,+),Ri%ein) du dv 
0 


n20 L,m=0 ¥. 


=, > (Rig), Ain(Rig,)”). 


On the other hand, the right-hand side becomes 
lim (Ry Qin ’ (R = RyRy ein) = lim (Pin ’ (Af = LD) gin) =; — 1. 


nro 


Hence, by equating both sides and dividing by ), , 


bam XO |S anralthey en een) 


Z,m=0 


090 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1967 
i p3 (Rip) OE (Rig), hiim = hjim) 


+ [ Cle) (dhinled Rig) Cu) du + (he)? , Hn Bt tea | (18) 
Thus, 


2 (1 _ 1) = = Pp Ajrimleistim(S; y Sx) 


=0 1,k=1 


ae PD Risia iva + hsim)(E;) ae [ Rytyn(Uju)Rin(U) du + tr (RsssvnBl a) | 
j=l 0 
< om, 
where (17) is used repeatedly. Hence,* 


Q0< I]. < om, 
+=0 


- (a) = (I rs) exp E 2d (1 7 1) ce) | , as., [Pi], 


where 
n(x) = liam. (2,RToin), [P.,P.] 7i¢=0,1,2,-°-. (19) 





Now, x°”(¢) has the following orthogonal expansion: 
e(t) = Lim. >> 1.(z)(Rig)?@, (Pil, OS Sr 
N00 4=0 
uniformly in ¢. Hence, there exists a subsequence of the partial sums 


7 (tz) (Rie;)‘” (t) which converges a.s. [P,] to x‘ (¢), uniformly 
in ¢t. Therefore, from (18) and (19), 


= li > (1 — 58) 
= Ojkimt (s0°”(s,) a > x(t; ae”, hjim + him) 
t,m=0 j 


ii. [ 2 (hina (u) du +(e, Hin o™) |, as. [Pi], 


which completes the proof of (22). 
* See Ref. 3, pp. 1653-1654. 
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Timing Recovery for Synchronous Binary 
Data Transmission 


By BURTON R.SALTZBERG 


(Manuscript received November 15, 1966) 


This paper analyzes different methods of adjusting the sampling time 
for detecting synchronous binary data, based on properties of the random 
data signal itself. The static error and the variance of the jitter of the 
resultant sampling instant are calculated where the effects of frequency 
offset, additive noise, signal overlap, and jitter of the reference source 
are included. 

The threshold crossing timing recovery system adjusts the sampling time 
an response to the tumes at which the data signal crosses the amplitude 
threshold. The samopled-derivative system uses the teme derivative of the 
signal at the sampling time to adjust sampling phase. It 1s shown that 
both systems lead to approximately the same amount of jitter tn the presence 
of noise and signal overlap for a given bandwidth of the control loop. 

An tmproved timing recovery system is presented which 1s constructed 
by adding correction signals to the sampled-derwative system. This system 
accounts for intersymbol interference in a manner that tends to set the 
sampling time at the point of maximum eye opening, where the error 
probability 1s minimum for the most adverse message sequence. 


I. INTRODUCTION AND SUMMARY OF RESULTS 


In synchronous polar binary data transmission, information 1s sent 
by serially transmitting either a basic signaling waveform or its negative 
at fixed time intervals. Modulation may be used to better fit the signal 
to the channel. At the receiving end, the signal is demodulated and 
filtered. The resultant baseband signal is sampled periodically, and 
the polarities at the sampling instants determine the output data. 
The choice of sampling time is critical for minimizing the error prob- 
ability due to intersymbol interference and noise, particularly when 
the signal has been subjected to sharp cutoff filtermg. The sampling 
time is best set by using some properties of the data signal itself. 
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The problem of timing is particularly acute in pulse code modulation 
(PCM) systems, where the accumulation of jitter in a long chain of 
regenerators frequently limits the allowable length of such systems. 
I’or this reason, previous studies of timing recovery have concentrated 
on PCM applications.’“* The use of a tuned circuit as the memory 
element is generally assumed, since this is commonly employed in 
PCM repeaters. 

This paper will concentrate on timing recovery for data transmission 
applications. The effects of multiple regeneration will not be considered. 
The recovery of timing will be accomplished by a feedback control 
system, such as a phase-locked loop. Different methods of generating 
the error signal for the control loop will be compared. 

The received signal, after demodulation and filtering, is of the form 


s(t) = p» af(t — BT — kT) + nd) (1) 
where {a,} is a set of independent random variables, each equal to 
+1 or —1 with equal probability. This may be assured by the use 
of a scrambler if the data source is itself not random. The basic signal- 
ing waveform is f(¢). The abscissa of f(t) will be adjusted for each system 
to be studied so that the desired sampling time of f(t) is t = 0. The 
quantity 8 is an unknown fractional time delay. Since we are not 
concerned with absolute time delay between transmitter and receiver, 
we will assume | 8 | S }. The additive noise is n(¢). 

The sampling wave which determines the times at which s(¢) is 
sampled may be represented by 


co 


qt) = , 6(¢ — nT — yT), (2) 


n=— 


where y is a phase that 1s generally time varying. 
The output data is determined by 


d, = sgn s(nT + yT), (3) 
where sgenv = v/|v |. Then 
d, = sgn [af@T — BT) + dS aufkT + yT — BT) + n(d). (4) 
k0 : 
For simplicity, the argument of the noise term is not made explicit 


since it is of no consequence. Assuming that f(7vZ — 87) is positive, 
then d, will agree with a, provided that 


— anf 2 An—J(kKT + yT — BT) + nd] < fT — BT). (5) 
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It is readily seen that the probability that (5) holds depends strongly 
on y — @. For each timing recovery system, f(t) is defined so that the 
desired value of y — B Is zero. 

_ The principal part of a timing recovery system is a phase detector 
which examines s(é) and g(t) and attempts to generate an error signal 
proportional to 8 — y. It 1s not possible to determine @ exactly from 
s(t) since the signaling waveform and the noise are unknown. This 
paper will consider different methods of forming this estimate of B — y. 

Another essential component of the system is the reference source 
which is used to generate the sampling wave. Its phase or frequency is ad- 
justed by the error signal in order to form a sampling wave of the proper 
phase. The error signal may be filtered prior to its use in shifting the 
phase of the reference source. The reference source may be a local 
oscillator whose natural frequency is set as close as possible to the 
bit rate. The reference source may instead be derived from transmitted 
pilot tones, in which case its frequency is exact, but it might have 
phase jitter of its own due to channel noise. 

Section IJ describes and analyzes a timing recovery system which 
uses a threshold crossing phase detector. This detector generates an 
error signal each time the signal crosses zero. The amplitude of this 
signal is proportional to the difference between the time of occurrence 
of the zero crossing and the time of the nearest sampling pulse, dis- 
placed by half a bit period. This system tends to choose a sampling 
instant which is midway between the mean transition times. 

The sampled-derivative phase detector is discussed in Section III. 
This device generates an error signal during each bit interval which 
is proportional to the time derivative of the signal at the sampling time 
multiplied by the signal polarity at that time. The sampled-derivative 
timing recovery system attempts to set the sampling time to coincide 
with the peak of f(d). 

The analysis shows that the performance of these two systems is 
very similar for a given open loop gain function of the control system, 
G(w). Approximations are made based on the assumption that the 
phase error is small and that G(w) is a narrowband low-pass function 
compared with the bit rate. 

Nhe systems fail if G(O) is finite and the reference source does not 
agree exactly in natural frequency with the bit rate. If the reference 
source has the correct frequency and a phase @, then a static phase error 


ss 6 B 
"14+ GO) (6) 


results in the sampling wave. 
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Much better performance can be achieved if G(w) has a pole at the 
origin. In this case the system is insensitive to the phase of the reference 
source. In the presence of a frequency error, Af, the static error is 


a= af = to |... (7) 


In addition to any static error, the sampling time will also Jitter 
about its mean value. The variance of this phase Jitter is the sum of 
several components, each due to a different cause. Jitter is produced 
by jitter in the reference source, by additive noise and by signal overlap. 
In the case of the threshold crossing system, jitter is also introduced 
whenever there is a static error. 

If the reference source has a jitter whose power spectral density is 
S,(w), then the output fractional jitter will have a variance equal to 


_1 S,(w) 
-{ 7+ter® 8) 


This indicates that high-frequency noise components must be removed 
from the reference source prior to its use for timing recovery. 
The jitter produced by the additive noise is 


Be = ; 2[F,(0) a RAT) és (9) 
“ — Tf(—7/2) — f(7/2)P 
for the threshold crossing system. /?,(t) is the autocorrelation of the 
noise and w, is the noise bandwidth of the closed control loop. 











__ Gw@) _ 
i+ G6) de. uo 
For the sampled-derivative system, the noise leads to jitter variance 
Bf 5. P!"(0) 
On = yaa) ° (11) 


In typical data transmission systems, (9) and (11) are similar in 
magnitude, and not very sensitive to the shape of /(é) if the noise is 
similarly filtered. 

The jitter variance duc to signal overlup is of the form 


o% = Aw T -+- Az(w.T)’, (12) 


er) 

3 = | 2 

Ono =—- ~ &) 
7 Jo 


where 
2 


_G@) | a. (13) 


1 + G@) 








TIMING RECOVERY FOR DATA TRANSMISSION 597 


If w,7’ and w.7’ are comparable and much less than unity, then the 
first term is usually much larger than the sceons 
For the threshold crossing system, es 


A, = ee eee 2 i (tkT + 1/2) — f(kT — T/2)] 
(26k + 1/2) + f(—kT + 7/2) -— f(—kT — T/2)] (14) 
and 
yn eee eee ree (arr? 
a0 (=2/2) — 7-7/2) 


+2 > f(kT + T/2)f(kT — T/2) — pe kf(kT + 1/2) 


— (er = T/MI(-EP + 2/2) — KEP — 7/2} (15 
Tor the sampled-derivative system, 
1 = , ? , P : 
A, = TO) a, PETE (AT) + f'(—-kT)] (16) 
and 
Aa = sper De BY EDY(-E1). a7) 


In both cases, A, = 0 if f(¢) is an even function, so the timing recovery 
systems are very sensitive to asymmetry of the basic signaling wave- 
form. The Jitter variances are again comparable for both systems. 
As may be expected, the jitter increases considerably as the filter used 
to shape f(¢) is made sharper. 

There is an additional jitter component for the threshold crossing 
system whenever there is a static error. Its variance is given by 


o = Bw, T'. (18) 


An example is provided in Section IV. A typical data transmission 
system using a distorted signal is studied so as to illustrate the mag- 
nitudes of the above quantities and to indicate the narrowness of loop 
bandwidth required for satisfactory performance. 

In this example it is also seen that ncither the threshold crossing 
timing recovery system nor the sampled-derivative system chooses a 
mean sampling time which Is very near to the time at which the eye 
pattern has its maximum opening. 
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Section V describes an improved timing recovery system whose 
mean sampling time coincides with that of the maximum eye opening. 
This system is constructed by adding correction signals to the sampled- 
derivative system in order to account for the effects of intersymbol 
interference on the mean sampling time. 

Finally, an outline of some extensions and modifications of these 
timing recovery systems is presented in Section VI. 


II. THE THRESHOLD CROSSING SYSTEM 


Most timing recovery systems make use of the instants at which 
the data signal crosses the threshold to alter the phase of the sampling 
wave. A block diagram of a typical threshold crossing timing recovery 
system is shown in Fig. 1. 

The principal part in this system is the threshold crossing detector. 
This device generates an error pulse each time the signal crosses zero. 
The amplitude of the error pulse is proportional to the difference 
between the time of occurrence of the threshold crossing and the time 
of the nearest pulse of the displaced sampling wave. The displaced 
sampling wave Is 


a) = qt- 7/2) = di s¢—-nT —T/2-y7T), (19) 


n=—00 


where y is the phase measured in fractional signal periods. If the axis 
crossing following the mth sampling time is displaced by a,,7’, 


s(mT + 7/2 +anT) = 0, l|anty| <i, (20) 









LOW-PASS 
FILTER 
G,(w) 






THRESHOLD 
CROSSING 
DETECTOR 


s(t) 









PHASE 
SHIFTER 


REFERENCE 
SOURCE 










Fig. 1— Threshold crossing timing recovery system. 
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then the error signal is 
ef) = Dd) Ki@n — y) 5 — mT), (21) 


where JX, is a gain constant. The error pulse has been represented in 
(21) as an ideal impulse function. Since the error signal will be passed 
through a narrow low-pass filter, the response will be virtually identical 
to that when a more realistic pulse of the same area is used. Similarly, 
the effects of variation of the position of the pulse within the interval 
may be neglected, since the low-frequency components of the error 
signal are substantially uneffected. 

We will now determine the threshold crossing time a,, as a function 
of the signal overlap and noise. Substituting (1) into (20) yields 


py AOn—-f(kT + T/2 — BT + a,T) +n) = 0. (22) 

If 
(2/2) + (-T/2) > De [IGP + 7/2) |+nQ (23) 
then a crossing will occur following the mth bit, if and only ifa,, = —@y41. 


When a,, = —Q@,+, and (23) holds, (22) may be written as 
Anlf(2/2 — BT + oT) — f(—T/2 — BT + anP)] 
+ 24 Om (kT + T/2 — BT +oanT) +n(t) =0. (24) 
If a, — 8 is small, we may approximate (24) by the first terms of its 
Taylor series expansion. | 
Omlf(T/2) — f(—T/2)] + amlom — BTA (L/2) — {(-T/2)] 
us m2 Om f(kT + 1/2) + n(t) 0. (25) 


Let the abscissa of the function f(¢) be adjusted so that 


f(T/2) = f(—-T/2). (26) 
Then define 


b= f'(—T/2) — f'(2/2). (27) 


We may now solve (25) for ap . 


Om BATE Lo Onl RE + 1/2) + n(Q)]. (28) 
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If we let 


d, cs ‘ if an = On+1 (29) 
iP if GQ, = ~An+1 


then the error signal (21) becomes 


et) = K, > a6 —-y+ ei 2 Ani (kT + T'/2) 


n=—0 


+ n(o)} 6(¢ — nT). (30) 


The error signal is passed through the filter G,(w) and then shifts 
the phase of a reference source. The reference source may be a local 
oscillator whose frequency is tuned as closely as possible to the signaling 
rate. Alternatively, the reference source may be derived from pilot 
tones which are transmitted along with the data. In the latter case, 
there is no error in the average frequency of the reference source, but 
its phase may be poorly related to that of the data signal and may 
also be perturbed by noise. In either case, the reference source generates 
a signal of the form 


rit) = >> b(t — nT — rT). (31) 
When a local oscillator is used, 
7r Aft + 6, AfT <1 (32) 


where Af is the frequency offset of the oscillator and 6 is an arbitrary 
constant. When the reference source is derived from pilot tones, 


r= 6+ n(2) (33) 


where y(t) is a zero mean random variable. 

The sampling wave is formed by shifting the phase of the reference 
source by an amount proportional to the value of the filtered error 
signal. 


q(t) = 3 6(¢ —nT — rT — KT), (34) 


n=—00 


where »v is the filtered version of e(t) and K, is the proportionality con- 
stant. Comparing (34) with (2), 


y=74+ Kv =7+ Kp g(t — x)e(x) dx. (35) 
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In the phase-locked loop type of control system, the frequency of 
the reference source is adjusted rather than its phase. This, however, 
is completely equivalent to integrating the error signal prior to phase 
adjustment. The block diagram is therefore valid for the phase-locked 
loop provided that the filter includes a pole at the origin. As will be 
shown, this pole is highly desirable and, in some cases, absolutely 
essential. 

Substituting (35) into (30) 


efi) = K, > 148 —r-—k, | g(t — x)e(x) dx 


n=-—0 


oe [Do da af(hE + 1/2) + noon 6¢— nT). (88) 


Let 

e(t) = ell (37) 

2 K, 

KK 

g(t) = = ll (38) 

e(t) = >> e(n) 6(t — nT) (39) 
and normalize the time variable so that 7’ = 1. Then (86) can be 
written as 


e.(n) = 2a.\6 —-T- > g2(n — k)e2(k) 


+ FLD aeaf( + 1/2) + no (40) 


A model of the threshold crossing timing recovery system which con- 
forms with (40) is shown in Fig. 2. This is not a time-invariant linear 
system because of the presence of the multiplier. 





a+) S, jan-kF (k +2) +n (t) 
. > 


Fig. 2— Model of threshold crossing system. 
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In Appendix A, it is shown that 
e(n) = B — + — g(0)e2(n) — Zs gn — kex(k). (41) 


This system may be readily analyzed either by means of the z-trans- 
form, or equivalently, by discrete Fourier analysis.” The discrete 
Fourier transform is given by 


io) 


X(w) = 2 a(n) exp (—jwn). (42) 


If G2(w) is bandlimited to | w | < a, which is approximately true in 
all cases of interest, then these transforms will coincide with the true 
Fourier transforms. 

The solution of (41) in terms of Fourier transforms is 


where 
9:(0) = = 1 ; Cra (44) 
The static error can now be found from 
1@) — B@) = G.@)E.) + 7@) — B&) (45) 
7) — B@) = 7 FLO __ fw) -6@l. 48) 


1 + g2(0) + G.@) 
In particular, if 8(¢) is a constant, 6) , and 7(¢) is given by (82), then 


(es) — By 8(0) = 2x GT Li Al OW) + (0 B)@)]_ 7 
10) — By = EO tage + 0 — 6 


~ 1+ g2(0) + G20) 


a oes A 1 + g2(0) 
~ J Aj dw |; ++ g2(0) + |... (48) 


If G,(O) is finite, then the first term is a steadily increasing error 
and the system fails. If Af = 0, the system does not fail, but a static 
error will be present due to the arbitrary phase values, 8) and @. It is, 
therefore, highly desirable that G(w) have a pole at the origin. In the 
presence of frequency offset, this pole is essential. In this case, only 
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the last term in (48) is nonzero, so that a frequency offset leads to a 
static error. The system is completely insensitive to the values of the 
arbitrary phases, 8) and 6, except during initial start-up. 
For the system to be stable, it is required that 1 + g.(0) + G2(w) 
have no zeros in the half plane Jm(w) < 0. In most cases, g2(0) « 1, 
so that the usual criteria for stability apply to G.(a). 


We now wish to calculate the variance of y, or the mean square 
jitter. From (40) and (41), and assuming the static error to be constant, 


é;(n) = e.(n) — @, = 2a, é -|- go(O)es — yy go(n — k)e;(k) 


— 7) + (9 | + 2(n) — Co (49) 


where 


a(n) =F dady So aynaflle + 1/2). (50) 
Let 
x(n) = 2a, 6 + 92(0)¢2 + a n(p + en) — @ (51) 
Then (49) may be written as 
ea(t) = a(n) — 2din() — 2d, D> guln = Wel). (62) 


The zero mean component of the output phase error is 


rin) = 00) = 7 = al) + DO gale — Heal (58) 
v0) =) + DO gor — Blo) — Aden). 4) 


The autocorrelation of (54) is 


DL LEI Pdig — b) + dul2digh + m — 1 + drench} 
= Rm) + Le Liga — Bgn + m—DRAk- 1, (55) 


where '(v) denotes the expected value of v and R,(m) = E[v(n)v(n + m)] 
is the autocorrelation of v. 
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In almost all cases of practical interest, g.(0) « 1. Then we may 
approximate 


a (56) 
and 

—_ f-¥ ~,_ x 

@=7T4 910) Bo — ¥. (58) 


Subject to these approximations, we can evaluate the discrete Fourier 
transform of (55). After some algebraic manipulation, the result ob- 
tained is 


[1+ Gx) 8,6) + |G) PRO = Se) + | Ge) FS), 69) 
where 
S.(4) = Do Ralon) exp (—jme) (60) 


is the power spectral density of v. 
The variance of y is calculated as 














ana, ne (61) 
G4) | S,@) 
du + 3 d 
; iGO Sl) FG) PO gay 
: 1 G Ee : 
Since G.(w) is narrowband, we may assume that 
a __ Gx(@) __ 
[ aan “de «1 (63) 








and therefore, 


o, 


GL) | 1 f° S,(w) 
mae “SiGe += | Teo ee (64) 
The second term indicates that low-frequency components of the 
reference source noise are attenuated while high-frequency components 
are not. Therefore, if the reference source is derived from transmitted 
pilot tones, it should be filtered to a narrow bandwidth before being 
used in the timing recovery system. 
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In Appendix B, S,(w) is evaluated for | w | « 1. This is the only region 
of interest when G.(w) is a sufficiently narrow low-pass filter. 


Sl) = +7 1RO) — R)) + Art Aw, 65) 
where 
Ar = js DS Ue + 1/2) - fb - 1/2) 
‘[2f(b + 1/2) + f(-k + 1/2) — f(-k- 1/2] (66) 


and 


As = sls} —4F/2) +2. 10+ L/h — 1/2 


— SF FUT 2) — 1-1/1 + 1/2) —(-k- 1/9} 6 








If we let 
_1 [*|_G)_ |’ 
, = =i EERE du) (68) 
and 
S> i 1 nae G',(w) : 








then (64) becomes 


= 2 
a, - C201 + ig [P,,.(0) = FC) Jo, + Aa, a Awe 


al SS, (w) = 
a iscar © 


This equation is given in the summary with the normalization 7 = 1 
removed. An application to a typical data transmission system is given 
in Section IV. 

It should be noted from (69) that w. will be unbounded unless G2(w) 
has at least two more poles than zeros. Good design of G.(w) requires 
that the second pole (assuming no zeros) occur somewhere in the 
vicinity of gain crossover. In this case, w. is approximately equal to w, . 

The first term of (70) indicates that the standard deviation of the 
jitter caused by frequency offset will be much less than the mean 
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of that error. The second term is the jitter produced by the additive 
noise. In typical data transmission systems, the signal is filtered such 
that &,(1) + 0. In that case, the variance of the jitter is directly pro- 
portional to the noise power. Of particular interest is the jitter produced 
by signal overlap. It can be seen from (66) that A, = 0 if f(¢) is an even 
function. In that case, the jitter variance is proportional to the cube 
of the system bandwidth, and can be made quite small by the use 
of narrow filtering. If f(¢) is not an even function, then the third term 
will usually be much larger than the fourth term. In either case, the 
jitter will greatly increase as the filter used to shape f(t) is made sharper. 


Ill. THE SAMPLED—DERIVATIVE DETECTOR 


An alternative method of adjusting the phase of the sampling wave 
makes use of the time derivative of the signal at the sampling times. 
Implementation of a sampled-derivative timing recovery system is 
about equally complex as a threshold crossing system. 

Iixcept for the manner in which the error signals are generated, 
the control loop is the same for both systems. Fig. 1 may be used to 
describe the sampled-derivative system if the delay in the feedback 
path is eliminated and a sampled-derivative detector is substituted 
for the threshold crossing detector. 

The sampled-derivative detector generates an error pulse during each 
bit interval whose amplitude is proportional to the time derivative 
of the data signal at the sampling time, multiplied by the polarity of 
the signal at that time 


e(t) = K; >> sgn [s(nT + yT)]s'(nT + yT) 6(t — nT). (71) 
However, the output data is generated by setting 
d, = sgn [s@T + yT)], (72) 


where 4, 1s the receiver decision on the nth bit. If the error rate is low, 


A 


Gd, = a, with high probability, and (71) may be approximated by 
e(t) & Kz >< a,8'(nT + yT) 6(¢ — nT) (73) 
if the effect of errors is neglected. 
Using (1), 
e(t) = K; >> anl Ant’ (kT + yT — BT) + n'()] 8(¢ — nT). (74) 
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The abscissa of f(é) in this case is adjusted such that the origin coin- 
cides with the peak of f(¢). 


f’(0) = 0. (75) 


If the phase error y — 8 is small, we may approximate (74) by the 
first terms of its Taylor series expansion 


e(t) YK; D> a,lany — 8)TF’’(0) + Dd An—if’ (KT) +n'()] (¢—nT). (76) 
As in the previous case, 

e(t) = >> e(n) d(t — nT) (77) 

y = K, a gi(t — aje(x) dx + 37 (78) 


and we normalize the time variable by setting 7 = 1. Equation (76) 
may now be written as 


e(n) = K,[(r — 8)f"O) + a, Dd, On~0f"(kK) + ann’ (2) 


+ K.f’’(0) = gin — et | (79) 


Let 
e,(t) = Se) (80) 
and 
g(t) = —K.K3f’’(0)9.(0). (81) 
Then (79) becomes 
elt) =B—7— yh) — D gn Dea, 2) 
where 
Qn / / 
yn) = 77) [2d dn—if’'(k) + n/(0). (83) 


Unlike the threshold crossing system, the sampled-derivative system 
is a time-invariant linear one when the phase error is small. A model 
of the system conforming with (82) is shown in Fig. 3. This model 
may be readily analyzed because of its linearity. 
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Fig. 3— Model of sampled-derivative system. 


yin) = D gn — Hest) + + (84) 
vin) = D> goln — 118 = yt) — 8] + +. (85) 


The mean error may be found from 


yw) [1 + Gs@)] = Gs@B@) + t@) (86) 
ls) — Bl) = -Ey@) = FEO. ($7) 


Equation (87) is identical to (46) if G3(w) is substituted for G.(w)/[1 -+ 
g.(0)]. All the comments of Section II concerning the static error and 
the desirability of G(w) having a pole at the origin therefore, also apply 
to the samplcd-derivative system. 

The variance of y when the mean error is constant will now be found. 


7102) = yin) =F = =D) gale — HlyH + 1 +09, 88) 


where 7(i) is again the reference source jitter. In terms of the power 
spectral densities, 


_|_ Gl) _ Sie) 
Sy) = | 1+ G() Sy) + 1 + Ga) [? (89) 
In Appendix C it is shown that, for | w |< 1, 

Swe eA ae (90) 


~ f’?0) 20) 


where 


A; = 





may HOU + (-] Qn 


=—0 
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and 
1 = Zr te 
A, = 7 2L, k’f'(b)f'(—k). (92) 
The variance of y is then 
2 1 f* 
c, = =a S,(w) dw (93) 
2 RYO) pp Lp. SO 
Cy = f’*(0) Wy + Aza, + A 4W3 + al | 1 ac G(w) ? dw, (94) 
where 
= G3) i 
Wy [Co do (95) 
and 
a. 1 __G;@) | ‘ 
He Ne ~[ 3 1 + G3(w) oe we 








There is a very strong similarity between (94) and (70). The last 
terms are identical, so that it is Just as important in the sampled- 
derivative system as in the threshold crossing system that high-fre- 
quency nolse components be removed from the referenced source prior 
to use for timing recovery. 

The jitter due to additive noise is proportional to the power of the 
derivative of the noise. The example in the next section illustrates 
that this is not serious if the noise is bandlimited to the same frequency 
range as the signal. However, if any high-frequency noise is allowed to 
enter the receiver beyond the signal filter, the jitter will be greatly 
increased. 

The jitter due to signal overlap is very similar to that of the threshold 
crossing system. If f(¢) is an even function, then its derivative will 
be odd, and A; = 0. Both A; and A, increase markedly as the spectrum 
of f(t) is made sharper. 

Tinally, unlike the threshold crossing system, there is no additional 
jitter term due to static phase error. This jitter component is eliminated 
because there is an error pulse generated during each bit interval. 


IV. AN EXAMPLE 


In order to illustrate the results of the previous two sections, the 
output jitter of both a threshold crossing timing recovery system and 
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a sampled-derivative system will be calculated for a typical received 
data transmission signal. 

The signal to be considered is one using half raised-cosine amplitude 
shaping which is distorted by linear delay distortion. For simplicity, 
it is normalized so that 7’ = 1 and the undistorted peak of the signal is 1. 





F@) = A@) exp [je@)] (97) 
1, 0O<w< 5 
A(u) = ; (98) 
pee 
2 2 ul om 
cos 5? 5 <a< 9 
o 2 
gw) = n°? wo > 0. (99) 


The outline of the ‘‘eye pattern” for this signal is shown in Fig. 4, 
along with the central portion of f(¢). The eye pattern is formed by 
superimposing the signals of all possible message sequences. Closing 
of the eye is due to signal overlap. 

The time of maximum eye opening is the optimum sampling time in a 
minimax sense. When such a sampling instant is chosen, then the error 


SDTC IM 


Fig. 4—Eye pattern outline and mean sampling times for a typical distorted 
data signal. 


TIMING RECOVERY FOR DATA TRANSMISSION 611 


probability in the presence of additive noise is minimized for the most 
adverse message sequence, where all adjacent bits overlap in a manner 
that subtracts from the bit being detected. 

The mean sampling times are also shown in Fig. 4 for the case of 
no static error. The mean sampling time of threshold crossing system, 
TC, is such that f(7C + 4) = f(T7C — 4). The mean sampling time 
of the sampled-derivative system, SD, coincides with the peak of f(d). 
It is seen that the threshold crossing system chooses a better mean 
sampling time than the sampled-derivative system, yet both systems 
miss the maximum of the eye opening by a Jarge amount. 

In order to compare the jitter due to noise, a particular noise spec- 
trum must be considered. Here it will be assumed that the noise 1s 
white noise which has been passed through a receive filter matched 
to the undistorted signal. In this case, the noise power spectral density 
is of the form 


Nw) = tA) (100) 
= Alcs) theo 
so that 
R,(0) = o; , (101) 
R,(1) = 0, (102) 
and | 
| ° w Aw) dw 
RY(0) = —~————- a. (103) 
i Ads 


From (70), the rms jitter due to noise for the threshold crossing 
system is calculated to be 


oy = 0.6150, Var . (104) 
For the sampled-derivative system, this quantity is calculated from (94). 
oy = 0.6120, Var . (105) 


The results are virtually identical. In either case, if c, = 0.1 (signal- 
to-noise ratio of 20 dB) and w, = 0.01, then the rms jitter due to noise 
alone will be 0.61 percent. It should be mentioned that several other 
signal pulse shapes were examined, and it was found that the Jitter 
due to noise was not very sensitive to pulse shape. 
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In order to calculate the jitter due to signal overlap, A, and A; 
were computed from (66) and (91). A. and A, were small enough to 
have negligible effect on the jitter. The resultant rms jitter is 0.287 Vw 
for the threshold crossing system and 0.265~/w, for the sampled- 
derivative system. If w, = 0.01, the jitter is 2.87 percent and 2.65 
percent, respectively. This is by no means negligible, and illustrates 
the need for very narrow filtering in the timing recovery control loop. 
Again there is little difference in the performance of the two systems. 

To observe the effects of asymmetry of the signal pulse, let us con- 
sider the same signal without phase distortion. Both timing recovery 
systems will then set a mean sampling time at the best point. The 
jitter due to signal overlap is greatly reduced since A, and Az; are zero. 
The computed values of A, and A, are 0.11 and 0.32, respectively. 
If w; = 0.01, then the rms jitter due to signal overlap is only 0.01 
percent for the threshold crossing detector and 0.03 percent for the 
sampled-derivative detector. Both values are completely negligible. It 
may be concluded from this calculation that both timing recovery 
systems are very sensitive to asymmetry of the signal waveform, both 
in terms of choosing the average sampling time and the resultant jitter 
about that time. 


V. AN IMPROVED TIMING RECOVERY SYSTEM 


It was seen in the previous example that both the threshold crossing 
timing recovery system and the sampled-derivative system led to 
average sampling times which differed considerably from the time of 
maximum eye opening. However, it is possible to modify the sampled- 
derivative system so that it does seek the time of maximum eye opening 
as the average sampling time. 

At any time ft, , the signal amplitude for the worst message sequence, 
assuming the current bit is 1, is 


D(to) = f(t.) — DB | f(to + KP) |. (106) 


In the region where the eye is open, D(t,) > 0, and the eye opening 
is equal to 2D(t,). If a sampling time ¢, is chosen such that D(t,) < 0, 
then errors will occur for some sequences even in the absence of noise. 

An experimental examination of the eye patterns of a large number 
of actual data transmission systems indicates that D(t)) is almost 
always a concave function of t, . Therefore, if ¢) is adjusted according 
to the gradient of D(t)), then the maximum of D will be found. 

It is therefore desired to generate an error signal whose average 
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value is 
c(t) = KD'(t — BT + yT). (107) 
If we again normalize the system so that 7 = 1, 
e() = Kif'@ — B+ 
—~ DMt+k—Bt+ysenft+k—-B+y]. (108) 


k¥0 

Equation (108) exists and is continuous except at those points ¢, 
where f(t, + k —B+y) = 0. 

Tig. 5 is a block diagram of a system which generates an error signal 
whose average is given by (108). The first term of (108) is the average 
error signal of the sampled-derivative detector discussed in Section III. 
The improved timing recovery system therefore, will consist of a 
sampled-derivative system with added correction signals. Enough cor- 
rection terms are used to account for those adjacent bits which may 
be expected to overlap significantly into the bit interval under con- 
sideration. 


SHIFT REGISTER 





s(t) SAMPLE 
AND HOLD 


DIFFER- SAMPLE 
ENTIATE He AND HOLD 


Z SAMPLE 

AND HOLD 
PHASE 
SHIFTER 

REFERENCE 
SOURCE 


Fig. 5—Improved timing recovery system. 
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The kth correction signal must have an average value 
e(tr) = Kf'(e + t) sgn fk + 4), (109) 
where 
ht pay: (110) 


This correction signal is formed by first generating an auxiliary 
signal x,(n) whose polarity is expected to agree with that of f(k + ¢,). 
The derivative of the signal at the current sampling time is multiplied 
by the polarity of the signal at a time displaced from the current time 
by k bit intervals. The result is multiplied by the polarity of the auxiliary 
signal to form the correction signal. 


é.(n + t,) = Ks’(n + 4) sgn s(n — k + 4) sgn a,(n + 4). (111) 


In order to account for overlap into leading pulses as well as lagging 
pulses, a fixed delay must be built into the system, as indicated in 
Fig. 5. This delay is equal to half the shift register length, so that the 
central cell of the shift register stores the polarity of the current bit, 
while the other cells store the polarities of preceding and succeeding bits. 

The auxiliary signal is formed by multiplying the value of the signal 
at the current sampling time by the polarity of the signal which pre- 
ceded this signal by & bit intervals. If k is negative, the polarity of a 
succeeding bit is used. The resultant is filtered by a narrow filter, 
H(w), to form the auxiliary signal x, . 


a(n) = >> h(n — m)s(m) sgn s(m — h), | (112) 


where the time displacement ¢, is ignored. 

If we assume that the error rate is low, as was done in Section III, 
then we may approximate sgn s(n) = a, with little loss of accuracy. 
Using this approximation and substituting (1) into (112), 


a(n) = 2) hin — m)[anm(t) + $(K) + One 2) Om-vf@)]- (118) 
The mean of x; is 
t, = f(k)H(0). (114) 
In Appendix D it is shown that the variance of x, is 


pee = [ | H(e) P E. - “oa | FQ) 7 du + Pyles) | dw, (115) 


where P;(w) is the Fourier transform of 


p(t) = fb — Of +d. (116). 
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If H(w) is sufficiently narrowband, then it can be seen from (115) 
that the variance of x, will be small. Then, if f(%) is not too small, 
it may be expected that the polarity of x, will agree with that of f(k) 
with high probability. 

We will now examine the mean value of the correction signal. Using 
the assumption of low error probability and substituting (1) into (111), 
the correction signal is 


e(n + t) = Kan >> dn-if'(n + bh) sgn a(n + hh). (117) 


In order to find the mean of (117), we make the approximation that 
x,(n) is independent of a,_, and a,_,. As a justification of this ap- 
proximation, note from (113) that 


Bp | Ane» One = F(A) + fRO) + fk — YAU — kb) (118) 
aCe (119) 
since h(n) «< H(0) for a narrowband filter. Let 
P, = Prob [sgn x, = sgn f(k)]. (120) 
Then 
e, & K(2P — 1)f'(k + t) sen f(k + 4). (121) 


When the magnitude of f(k + ¢,) is sufficiently large, P ~ 1 and the 
mean of the kth correction signal is approximately the desired value. 
In the vicinity of a zero of f(k + ¢,), 4 < P < 1, and the correction 
signal will at least have the correct polarity, although not the correct 
magnitude. Ideally, the correction term should be a discontinuous 
function of ¢t, at a zero of f(k + ¢,). The actual correction signal will 
have a mean value which is continuous, but the sharpness of change in 
the vicinity of a zero will increase as H(w) is made narrower. 

The rms jitter of this timing recovery system is extremely difficult 
to evaluate because of the presence of many nonlinear operations. 
However, this jitter may be expected to be much greater than that 
of a sampled-derivative system, since each correction signal may be 
expected to introduce jitter of the same order of magnitude as the 
main error signal. Narrow filtering in the control loop is therefore 
essential. 

The mean sampling for the example of Section IV when this improved 
timing recovery system is used is shown in Fig. 4 as “J//’’. It is seen 
that the time of maximum eye opening has been found. For this example, 
only one leading and one lagging correction term were sufficient to 
choose this mean sampling time. 
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VI. EXTENSIONS 


Kach of the timing recovery systems described here may be modified 
to work with m-level digital data signals instead of only binary formats. 
Since in multilevel systems the eye is much narrower, the effects of 
static timing error and jitter are much more serious. 

A threshold crossing detector may be constructed which generates 
an error pulse whenever the signal crosses any one of the m — 1 thresh- 
olds. During any signal interval, any number of error pulses between 
zero and m — 1 may be present. Such a system is extremely difficult 
to analyze, but has been found to work well in practice, provided 
that some auxiliary means is used to correct the mean sampling time. 

The sampled-derivative detector is very easily extended to multi- 
level systems if the signal derivative is multiplied by the output symbol 
value in forming the error signal. If the a,’s are scaled so as to form 
a set of unit variance, then the analysis of Section III applies directly. 

The improved timing recovery system is modified for use with multi- 
level signaling in a manner similar to that of the sampled-derivative 
system. However, the signal margin against noise for the worst message 
sequence is now 


D(t) = f(to) — (m — 1) > | f(t) + kT) |. (122) 


Comparing this criterion with that of (106), it is seen that each of the 
correction signals must be weighted by the quantity m — 1 in order 
to find the time of maximum eye opening. 

IXxtension of these techniques to partial response systems is also 
straightforward. Since the modifications depend on the particular 
partial response system used, a description will not be presented here. 

All of the systems analyzed here used linear control loops. This 
permitted the calculation of jitter variance in terms of loop bandwidth. 
However, the implementation of these systems may frequently be 
simplified considerably by using nonlinear control systems. A particular 
method which has met with practical success uses the polarity of each 
error pulse to adjust the sampling phase by a small fixed increment. 


APPENDIX A 


Evaluation of e,* 
Equation (40) is of the form 
€2(n) = 2a, ¢ = > go(n = eat) | ’ (123) 


—cs 


* The approach used here was suggested by J. E. Mazo. 
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where d, = 0 or 1 with equal probability and are independent. The 
random variables c, are uncorrelated with the d,’s. 
From (123) and the properties of d,, , 


€o(n) = d,e2(n) (124) 
and 
d,e.(k) = ge.(k), if n> k. (125) 


We first find the average of (123) over all c, and all d,, k < n. This 
average of a random variable v will be denoted by (v), while the average 
over all c, and d, will be shown as 3. 


(ealn)) = die, — 2 DO gum — bl d,ex(X)). (126) 


Using (124) 
n—-1 


(é2(n)) = 2d,e, — 29(0)(€2(n)) — 22 gain — k){d,e2(k)). (127) 


The only random variable in (127) is d, . The overall average, e,(n), 
is therefore the average of (127) over d, . Using (125), we obtain 


a) = 6 — 294(0)ea0) — Salm — Hex (128) 
alt) = & — ga(O)eat) — DO gale — Hex (129) 


which 1s the result shown in (41). 


APPENDIX B 


Evaluation of S,(w) 


From the definition of d, in (29), we may express 2d, as 


2d. = | — 6,044. (1380) 


Then (51) may be rewritten as 
= ] 
(MN) RI Any srl2 + 2M) HF (An — Aner)N(), (131) 


where the term g2(0)e, has been neglected. 

We wish to evaluate the power spectral density of x. It will first be 
shown that the approximation of x given in (181) is zero-mean. Since 
the a,’s are zero-mean and independent, and the noise n(¢) is zero- 


618 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1967 


mean and independent of all other random processes, then the first 
and last terms of (131) are zero-mean. The mean of z is readily found 
after substituting (130) into (50). 


a(n) = 5 (tn — aes) Dr deaf + 9) (132) 


z(n) = 0. (133) 


The mean value of x(n) as given in (131) is therefore zero. It may 
similarly be shown that the three terms of (131) are mutually uncor- 
related. The autocorrelation of x is then 


R,(m) ma QnOn+14n+m4n+m+1 e> + R(m) 


+ GH Gaon Gem Rm) — 034) 


R,(0) = & + RO) + SRO) (135) 
R,(e1) = R.(E1) — 7 Ry(cel) (136) 
Rm) =R{(m, m0, #1. (137) 


From (182), 


1 8 ee ee 
R,(m) = b? > > (a, a Cina) Gaga - Oitnikt) Oncaea 


k0,-1 1+0,—1 
fk + a)fl + 2) (188) 
2 


RO) =F Dy Pe +¥ (139) 
RI) = el D e+ MG - 9) + 1-VO! (140) 


R,(em) = 45 [fm +3) — Hm — Dlif(—m + 4) — f(—m — BI, 
m0, #1. (141) 


The power spectral density of x can now be calculated for | w | < 7z. 


2] 


S,(@) = p2 R,(m) exp (—jwm) (142) 


=— 00 


S.(w) = R,(0) + 2 > R,(m) coswm. (148) 
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Since x will be passed through a very narrow filter, we are interested 
in the spectrum of x only in the region | w | < 1. In this region, (148) 
may be approximated by 


S.(e) & R,(0) + 2 PB R.(m) — oo »? mR,(m). (144) 


This approximation is valid if the third derivative of S,(w) is bounded.” 
This will be true if R,(m) decreases as O(m7*). If the Fourier transform 
of f{(é) is continuous, then f(m) decreases as O(m™’). In this case, it 
can be seen from (137) and (141) that R,(m) decreases as O(m™°), 
so that the approximation (144) is valid for | w | <« 1. 


Sul) VE +S {2.00 — RQ) + D PE+D 
~ Y ik + dik — ) — (-DIG 


bol 

NS 

— 
—S 


+ Yom + 8) = fom = Dlif(—m + 9) = f(—m 

+o De + De — 9 + DIO 

— ¥ migim + — fm — Hli(—m — ) — fm — pih- gs 
After some manipulation, and using (26), (145) may be reduced to 
s.) = + 24{R@ — RQ) +2 D ge+H - 1-9 
12k +) + f(-k +8) — 1k — ah 


op {af +2 D0 fk + Hib — 4) 


=—0 


ae 


I 


— 3 wie +) — 1 — DI—k +) — HE DI} C49 


APPENDIX C 
Evaluation of S,(w) 
We wish to find the power spectral density of 


On , , 
y(n) = "Oy [ Ss nif’ (k) + n’(0)). (147) 
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The autocorrelation of y is 


none Tay [DO eypentin stromal COND 
+ G,Oninn’(i)n'(t + m)} (148) 
pas oe 72 — pr 
R,(0) = 777™(6) [df (k) — R7’(0)]. (149) 
For m # 0, 
ee 
R,(m) = 770) f’(m)f’(—m). (150) 


Under the same conditions stated in Appendix B, the power spectral 
density of y may be approximated in the region | w |< 1 by 


S,(w) © F,(0) + 2 > R,(m) — w > mR,(m). (151) 


Using (75) we obtain 


ia) 


Slo) & wag (RYO) +S POU® + 1-H) 


~ it S Kyaw} 089 


APPENDIX D 


Evaluation of o%, 
The variance of x, 1s the mean square value of the zero-mean random 
process 


x,(n) — == BS h(n ~ m)a,—,{n(t) + ps On-pf(p)]. (158) 


Since the noise and the message are independent, the variances 
of the two components of (153) will add to form the total variance. 
The variance of x, due to the noise 1s 


of = DY DU A(n — mh(n — Q)an—4Oq—1(tm)(E,) (154) 
a = do h?(n — m)n>(t) (155) 


a; = o, [- | H(w) |° dw. (156) 


| 
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The total variance of x, 1s then 
ork a a; a Dd ys Ds 2 h(n | m)h(n _ DO kAg—kAm—pAg-+ 
a pe¥k rk 
‘{@)f(M. (157) 


Nonzero contributions arise from those terms where gq = m and 
= p and from those terms where r = 2k — pandq = m+k — p. 


or = ort DY Dd hn — mf) 


m p#k 


+ hin — mh(in —- m — k + p)f(p)f(2k — p)). (158) 


Any one of the terms of (158) is small compared to the sum. We 
may therefore approximate (158) by including the missing p = k terms. 


owoitzs | | Hw) [dof |F@) Pdo+os, (159) 
where 


= » 2 h(m)h(m — k + p)f(p)f(2k — p) (160) 


and it is assumed that H(w) is bandlimited to | w| < 7. 


= aaa [fff mon@rory ¥ ex lime +a) 


- >> exp [jp(u + v — y)] exp [jkQy — w]dwdudvdy (161) 


Pi z i | | [ ; ON ONUOV LCE Cea Cee ee 
-exp [jk(2y — u)] dw du dv dy (162) 
i= Gs | H@H(-e)de [Ply — FQ) exp ky — a) dy (163) 


i= a] |H@) Pde | Fe — 4) exp [-ike — W)] 


-F(y) exp Gky) dy. (164) 


The second integral in (164) may be recognized as 27 times the 
Iourier transform of 


pit) = fk — Ok + b) (165) 
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so that 


gael : ; ay PUG) de. 266) 


2 De 


Substituting (166) and (156) into (159) 


o. = -- Z | Hw) |? E +- os [- | Fu) |? du + Put) | dw. (167) 
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Modulation of Laser Beams by Atmospheric 
Turbulence — Depth of Modulation* 


By M.SUBRAMANIAN and J. A. COLLINSON 
(Manuscript received September 26, 1966) 


We have studied the fluctuations produced in a laser beam by atmospheric 
turbulence over transmission paths up to 2400 feet long as a function of 
size of receiving aperture, range, and atmospheric conditions. The depth 
of modulation decreases rapidly with increasing size of receiving aperture 
for apertures smaller than the direct beam. It does not go to zero, however, 
but rather levels off at an approximately constant, finite value for apertures 
larger than the direct beam. 

When all of the direct beam is eoliedied the depth of modulation varies 
approximately with the 3 power of range from about 100 to 2400 feet, 
the largest range used. At ranges less than about 100 feet, however, the 
dependence 1s consistently much less than 3. These results are independent 
of weather conditions, of tume of day, of local conditions along the path, 
of whether the transmitter 1s inside or outside a building, of a twofold change 
an diameter of the launched beam, of whether the range 1s a straight pass 
or 1s multiply-folded, and of mirror separation in the multiply-folded 
arrangement. The 3 dependence is consistent with near-field scattering 
theory and leads to an estimate for the lower bound of the effective scale 
size of turbulence of & centimeters. 

The depth of modulation, however, arenas sensitively on atmospheric 
conditions; in a time of the order of seconds the value can change as much 
as an order of magnitude. We have systematically measured depth of 
modulation of the direct beam simultaneously with wind velocity and 
variability, temperature gradients, and time of day. No simple dependence 
on these variables was found. 


I. INTRODUCTION 


From the point of view of communications, one of the serious effects 
of the atmosphere on propagation of laser beams is the fluctuations In 


* Part of this paper was presented at the 1966 International Quantum Elec- 
tronics Conference, April 12-15, 1966, Phoenix, Arizona, 
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the received signal caused by variations in the dielectric constant of 
the air. An important measure of the fluctuations is their power spec- 
trum. Hogg! first measured these spectra and obtained an exponential 
distribution with a baseband width of the order of a few hundred 
cycles. Hogg used a multimode 6328 A laser and a range of 2.6 km. 
The receiver, 5 cm in diameter, was located in the center of the re- 
ceived beam, which was about 25 em in diameter. Hogg observed that 
an increase in the angular beamwidth of the source caused an increase 
in spectral width. 

Hinchman and Buck? measured the low-frequency fluctuations in a 
6328 A laser beam at distances of 9 and 90 miles. Their receiver, 3 
inches in diameter, collected a very small fraction of the total beam 
power. They observed a very large depth of modulation. The spectral 
density of the fluctuations was found to decrease with increasing fre- 
quency up to 50 Hz, the highest frequency measured. 

Subramanian and Collinson® propagated a single-mode, diffraction- 
limited 6328 A beam and examined the dependence of the spectrum 
on a variety of parameters. The transmitted beam diameter was 
changed from 1 to 38 mm, beam divergence was adjusted by focusing 
the telescope employed, ranges of 120 and 360 meters were used, and 
the receiver aperture was varied from much smaller to much larger 
than the received beam size. The spectrum, which had an exponential 
distribution (in agreement with Hogg), was independent of these 
variables within experimental error. The width of spectrum, however, 
was sensitive to atmospheric conditions. In general, the spectrum be- 
came wider as refractive gradients along the path became larger. For 
example, temperature gradients (caused by the sun) and pressure 
gradients (caused by turbulent wind) systematically gave broader 
spectra. The total width of the spectrum (above detector noise) varied 
over a total range of 60 to 1000 Hz, with typical value of a few 
hundred hertz. Thus, the spectral width was about the same as Hogg’s, 
although the distance was about an order of magnitude less. 

Buck‘ took additional data at more distances on the same 90-mile 
range, the smallest distance being 550 meters. Although the analytic 
shape of the power spectra was not given, it is nevertheless significant 
that the spectra shown all reached cutoff at about 200 Hz, and he 
stated that there was no systematic variation in the spectra when the 
detector aperture or path length was changed. Buck commented that, 
with a very large aperture, a noiseless de signal is obtained, but this 
was not found by Subramanian and Collinson.? Thus, the three sets of 
observers at three locations have found a characteristic width of 
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about a few hundred hertz, rather independent of the experimental 
arrangement, and, in particular, apparently insensitive to changes in 
distance from 120 meters to 145 km. 

Another important measure of the fluctuations in the received signal 
is depth of modulation, or the ratio of the rms power in the fluctua- 
tions to the average beam power. While spectral width may be inde- 
pendent of the experimental parameters, one expects the depth of 
modulation to show a strong dependence, especially on distance. The 
present work was undertaken to establish the nature of the dependence 
of modulation depth on such variables as distance, receiver aperture, 
and atmospheric conditions. 


II. EXPERIMENTAL ARRANGEMENT 


Since beam diameter will vary with distance (due to diffraction as 
well as atmospheric refraction) and since modulation depth presum- 
ably will depend on receiver aperture and on distance, the experiment 
must be arranged to allow adequate separation of these variables. For 
changes in distance to be meaningful, the receiver must bear some 
appropriate, uniform relation to the beam size regardless of distance. 
The simplest approach is to use apertures always larger than the 
direct beam, since it is known? that the fluctuations do not then 
clisappear. 

In order always to collect substantially all of the beam, the distance 
used should not be too large. Otherwise, the beam will be large and a 
receiver of an inconvenient size will be needed. With horizontal paths 
near the ground, it is commonly observed that atmospheric refraction 
produces angular spreading of the order of 10~* radian, so paths miles 
long would imply beams feet in diameter. Moreover, it is desirable to 
be able to change the distance essentially continuously, and ranges 
where this can be done beyond a few hundred feet are not easily 
obtained. 

Such short distances imply a very low level of atmospheric modula- 
tion with some values of modulation depth lower than 0.1 percent. 
This means, in turn, that the amplitude of the noise of the laser must 
not exceed about 0.01 percent. (In all cases, percent modulation is 
defined as 100 times the ratio of the rms power of the fluctuations to the 
average power). The laser used was designed’ for high intrinsic fre- 
quency stability, and, when properly operated, it has the necessary 
amplitude stability as well. The RF power supply must be well regu- 
lated, and de (rather than 60 Hz) power must be used on the filaments 
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of the power supplies. At the frequencies of interest, acoustically- 
coupled noise is substantial in the average laboratory and must be 
reduced. 

An easy and highly effective method of isolation is to seal the entire 
laser (at one atmosphere) in a gas-tight container. This was done with 
a glass bottled fitted with an optical-quality window at one end and a 
polished flange at the other end. The flange made an O-ring seal to a 
Lucite plate. (It is clear that pressure fluctuations cause considerable 
amplitude noise, since enclosures which do not form vacuum quality 
seals are not as effective.) It is helpful to place a felt pad between the 
laser and the bottle and Isomode (corrugated rubber) pads between 
the floor and the legs of the table on which the laser and bottle rest. 
With this arrangement, the amplitude modulation of the laser could 
be maintained at or below 0.005 percent. 

The laser used> was also single-mode and RF-excited in order to 
further ensure that measurements did not include any spurious noise. 
Amplitude fluctuations can appear in the output of multimode lasers 
as a result of mode competition and self-beating effects. Hodara has 
calculated® and measured’ the excess noise caused by mode-beating in 
multimode lasers. He observed that many of the reported discrepancies 
in measurements of laser noise probably arose because multimode 
lasers sometimes were used. 

While our method of mounting the transmitter is not massive, never- 
theless no detectable variation in beam pointing occurred during an 
experiment, and the arrangement has the advantage that it could be 
readily modified. As will be seen, the variety of necessary experiments 
required flexibility in both the transmitter and the receiver. The trans- 
mitted beam emerged from the room through a selected Lucite window. 
Ordinary plate glass caused noticeable refraction of the beam, but 
some parts of new panels of 4-inch thick Lucite produced no measur- 
able distortion of the beam. Many “A-B” experiments were made, and 
no difference could be found between modulation results with a Lucite 
pane and the results with an open sash. 

The range was located on the flat roof of a two story building at the 
Whippany location of Bell Telephone Laboratories. The roof surface 
was asphalt and gravel, but usually most of the path was over a 
wooden catwalk whose surface was 3 inches above the roof. The beam 
path was 3 feet above the roof and ran 30° east of north to a total 
available range of 330 feet. 

The receiver took a variety of forms, and it will be best to give each 
experimental arrangement with the corresponding results. However, 
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many of the details of a representative setup can be seen in Fig. 1. 
The equipment cabinet on casters gave the mobility which was re- 
quired for rapid change of propagation distance. Atop the cabinet is a 
section of triangular tower used when a long focal length, large aper- 
ture lens provided the receiving aperture. In the arrangement shown, 
a diffraction limited, 6-inch diameter, 8-foot F.L. lens was mounted in 
the right end of the tower section. The detector, an RCA 7265 photo- 
multiplier, was equipped with a 3 A wide-interference filter. The detec- 
tor housing appears in the left end of the tower section. 

Just to the left of the tower section is a bank of 6-inch diameter, 
4-wave flat mirrors supported by a bench which crosses the catwalk. 
Such mirrors were used to fold the beam and provide transmission 
paths up to 2400 feet long. The bench rested directly on the roof, and 
this provided adequate stability of the multiply-folded optical path. 

On the bench in the foreground can be seen an RCA vacuum tube 
voltmeter, used to measure the de voltage across the photomultiplier 
load, and a Hewlett-Packard 403A ac voltmeter, used to measure the 
rms ac voltage. The bandwidth of the 403A is 1 Hz to 1 MHz. To the 
best of our knowledge, this is the only ac meter that has such a low 
frequency response. Such response is important since the fluctuations 
are exponentially weighted toward the low end of the few-hundred 
hertz band. Since instantaneous voltage output from the photomulti- 





Fig. 1— Receiving station. 
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plier is proportional to the instantaneous incident optical power, it 
follows that percent modulation is given by 100 times the ratio of the 
ac voltage to the de voltage. 

Since the fluctuations do not disappear? when all of the direct beam 
is collected, it is necessary to show that possible sources of spurious 
noise are negligible. The laser, for example, was checked periodically 
to insure that transmitter fluctuations were no larger than about 0.005 
percent. Another possibility was that the fluctuations resulted from the 
product of the time-varying intensity profile of the beam and the 
uneven sensitivity profile of the photocathode. Also, the noise might 
have resulted from the dancing of a small beam (order of a millimeter 
diameter) over the sensitivity profile of the photomultiplier. Two 
experiments were made to assess these possibilities. 

First, a photovoltaic cell (Hoffman 110C) which was one cm square 
was interchanged in “A-B” fashion with the photomultiplier. When the 
cell was used, the one-cm beam was focused to about one millimeter in 
diameter, and all of it therefore was collected by the cell. The spatial 
variation of sensitivity of the cell was orders of magnitude smaller 
than that of the photomultiplier, and the measured fluctuations were 
the same as with the photomultiplier, within experimental error. 

In the second check, a comparison was made between the fluctua- 
tions obtained when a 1-cm diameter beam fell directly on the photo- 
multiplier, when a 4-inch diameter, 60-inch focal length diffraction 
limited lens focused the beam to about a 4-cm spot, and when the 
lens focused the beam to a spot about 4 mm in diameter. The fluctua- 
tions were the same in the former two arrangements. In the last case 
of the $-mm spot, however, when there were mechanical disturbances 
of the lens-photomultiplier assembly large enough that dancing of the 
spot on the detector was visible, the level of fluctuations was measur- 
ably higher. Since the photocathode had a sensitivity structure with a 
scale size of the order of a few millimeters, it seems clear that the 
excess noise was caused by the random scanning of the small spot over 
the spatially varying sensitivity of the detector. As a result, in all of 
the work which employed a collecting lens, the spot was deliberately 
defocussed to about 4+ to 1 cm in diameter, and mechanical disturb- 
ances of the receiver were avoided. 

The first and dominating result of any measurement of depth of 
modulation is that the value changes steadily with time. Another way 
of stating this is that the fluctuation spectrum extends well below 
1 Hz, the low-frequency limit of the ac meter, and the ac value changes 
steadily. The consequence is that this temporal change in modulation 
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is unavoidably combined with the dependence of modulation on the 
measured variables (such as distance) when the data are taken at 
different times. Thus, when distance is changed by moving the receiver, 
the value obtained changes because of variations both in distance 
and in time. 

The seriousness of this depended on the degree of temporal insta- 
bility of the modulation. Typically, the measured value of percent 
modulation varied by a factor of about two or perhaps three in a 
period of a few minutes. This was tolerable, but it meant that for a 
determination of distance or aperture dependence to be meaningful it 
was necessary to average the results for a large number of runs. For 
this reason it was important to arrange the experiment so that succes- 
sive readings could be taken quickly. In general, it was possible to 
move the receiver and make a reading in about two or three minutes. 
Thus, a run involving five points took about 10 minutes. On many 
occasions, however, the temporal instability of modulation was severe, 
and the value might change by an order of magnitude within a few 
minutes. This is comparable with the total change produced by the. 
variations in distance and aperture, and useful data could not then be 
obtained. Most of our results were taken at night, a few hours after 
the sun had gone down. During this period, changes in weather condi- 
tions were relatively small. Besides, the alignment of the receiver 
could be made very rapidly at night, hence a number of runs could be 
taken in quick succession. 

An alternative which would circumvent this problem would be to 
divide the beam with beam splitters into receivers at each distance or 
aperture value. Modulation then would be measured simultaneously at 
all the receivers. However, this would have required more extensive 
facilities than were readily available. 


III. RESULTS—APERTURE DEPENDENCE 


It should once again be emphasized that while the fluctuation spec- 
trum may be insensitive to size of recelving aperture, one expects the 
depth of modulation to depend rather critically on it. In particular, if 
the fluctuations are produced entirely by variations in power collected 
by a finite aperture, one might expect the depth of modulation to 
decrease as larger apertures are employed. With a large enough re- 
ceiver, the fluctuations should then go to zero. 

Depth of modulation was measured with a beam which appeared to 
the dark-adapted eye to be about 3 to # inch in diameter and with 
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Fig. 2— Aperture dependence: (a) schematic, (b) data. 


receiving apertures from } to 3 inches in diameter. Fig. 2(a) shows the 
schematic of the arrangement. Aperture dependence was thus measured 
at various distances. At each distance, several successive data runs 
were taken and averaged. The results are shown in Fig. 2(b). 

There are two curves for the 25-foot distance, one with considerably 
higher percentage modulation than the other. The former was taken 
before a rain storm and the latter immediately after the rain. Al- 
though the rain appeared to have a decided effect in this case, firm 
conclusions should not be drawn since this was a single observation. 
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(Nature did not present us with more than one such opportunity.) 
Dependence on the atmospheric conditions will be discussed more 
fully below. 


IV. RESULTS—-RANGE DEPENDENCE 


The dependence of depth of modulation on range was measured, 
using apertures large enough to collect all of the direct beam. The 
aperture dependence [See Fig. 2(b) ] with small apertures is sufficiently 
strong that a meaningful range dependence would be difficult, or im- 
possible, to obtain using apertures smaller than the beam. Generally, 
the aperture used was at least twice the size of the direct beam as it 
appeared at night (i.e., to a dark-adapted eye). Thus, all data for 
range dependence were obtained in the region in which the value is 
sensibly independent of aperture size. 

The schematic of the first experiment is the same as given in Fig. 
- 2(a) except that a 20-power telescope was used instead of the 9-power 
one. The averages of all the data are plotted in Fig.-3 on log-log 
coordinates. (The bars are averages of the mean deviations in the 
data for each night. The deviations are indicative of the unavoidable 
uncertainty caused by changing atmospheric conditions). 

The results from 100 to 300 feet suggest that depth of modulation 
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Fig. 3— Range dependence—single pass—(25-300 ft); 


—night—4 runs, 
5/20/65 night—8 runs, 
5/26/65 night—6 runs, 
5/28/65 night—6 runs, 
6/24/65 night—2 runs. 


632 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1967 


goes approximately with the 3/2 power of distance. The data at 25 
and 50 feet seem to indicate an end effect not intrinsic in the atmos- 
phere. It suggests that the source was noisy, producing spuriously 
large values at small ranges. To examine this possibility, the data 
were replotted on linear-linear coordinates, and the curve was extrap- 
olated to the y-intereept. This yielded an apparent zero-distance 
modulation of 0.024 percent. This is about five times the laser noise 
of 0.005 percent measured inside a quiet room, so that the laser can- 
not be the dominating cause of the change in curve shape at short 
distances. 

Another possible explanation is that local conditions affected the 
atmosphere differently along the beam. Such variations could have 
been caused by a row of four large exhaust blowers along a line about 
30 feet east of the range. On the night of June 24, 1965, we arranged 
to turn off all the blowers. Two runs were taken, and the results 
showed the same distance dependence as with the blowers on, so 
closely, in fact, that the data were simply included in the final curve 
of Fig. 3. 

Another possibility is that local conditions affected the atmosphere 
differently near the transmitter room than far away from it. For 
example, the room itself may well have changed the turbulence of 
the wind. We therefore set up the transmitter outside the room and 
30 feet away from the room (which is 8 by 14 feet). If local condi- 
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Fig. 4— Range dependence: single pass (10-200 ft) with transmitter located 
outside the building; 
8/3 /65—day—6 runs, 


8/3 /65—night—6 runs. 
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Fig. 5—- Range dependence: single pass (0-120 ft) with transmitter located outside 
the building and telescope changed to 9 power; 


7 /20/65—night—4 runs, 
7/21/65 day—4 runs, 
7/21/65 night—4 runs, 
7/22/65 day—d runs. 


tions were the cause, the entire curve should displace to the left on 
the abscissa by 30 feet. 

The schematic of this arrangement is identical to that shown in Fig. 
2(a) except that the laser is located outside the building with a 20- 
power telescope. The results are given in Fig. 4. The curve did not 
shift along the abscissa and shows generally the same behavior at 
the same distance from the transmitter. Whatever the cause of the 
near-distance behavior of the modulation, it seems to have moved 30 
feet with the transmitter. 

Because of the surprising behavior of the depth of modulation at 
short ranges, more measurements were made, still with the transmitter 
outside, but now adding some very short distances. The 20-power 
telescope was replaced with one of 9 power, giving a reduction in 
diameter of transmitted beam from about 2 to about 1 cm. This was 
done to examine what effect beam divergence might have on the dis- 
tance at which the knee of the curve appears. Data were taken during 
both day and night. The daytime data were found to be not signifi- 
cantly different from the night data. All values were averaged and are 
plotted in Fig. 5. It appears that the relatively large values of modu- 
lation depth which yield a small slope at short distances persist at 
distances as small as 24 feet from the transmitter. (The laser noise 
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level of 0.005 percent is measured at distances of this order, but in 
a laboratory where the air turbulence is low.) 

There appears in Fig. 5 to have been a significant change from 
the distance dependence displayed in Figs. 3 and 4. The knee of the 
curve seems to have moved to larger distances such that the 80-foot 
value now is aligned with the dependence at short distances. Since the 
change in beam size was an obvious potential explanation, and since 
the principal experimental difficulty still was the large temporal varia- 
tion in modulation, the following experiment was conducted. The 
laser was mounted outside and 30 feet away from the transmitter 
room, and the beam was split into two beams of approximately equal 
intensity. One beam was transmitted with the 20-power telescope, the 
other with the 9-power telescope. The beams were parallel and 20 
inches apart. At each distance, the modulation was measured on both 
beams before changing distance. Three runs were made on one night. 
This did not give enough data to define smooth curve shapes, but it was 
enough to show that the curve shapes were nearly identical for the 
two beams. The change in Fig. 5 from Figs. 3 and 4, therefore, cannot 
be attributed to the change in telescopes. 

Having explored the behavior of modulation at short distances, we 
now turned to distances larger than 300 feet and in particular to the 
question whether the 3/2 power dependence of modulation depth on 
distance would continue at large distances. Fig. 3 suggests that modu- 
lation depth at 300 feet may be a little lower than a 3/2 power de- 
pendence would imply. Results for three of the five nights summarized 
in Fig. 83 showed a rather pronounced reduction in the value expected 
at 300 feet by extrapolation from shorter distances. 

In order to obtain much longer paths in the available space, we now 
folded the beam back and forth over the range, as shown in Fig. 6(a). 
The laser again was mounted in the transmitter building, and the 20- 
power telescope was used in order to reduce beam spreading by diffrac- 
tion over these longer distances. At the maximum distance of 2100 
feet, the spot generally was about two to three inches in diameter. 
The mirrors used were always substantially larger than the beam. 
In Fig. 6(a), the first two mirrors were four inches square, and the 
last four were six inches in diameter. All the mirrors were flat to 
a quarter-wave and front-aluminized. Fig. 1 is a photograph of the re- 
ceiving end of this arrangement. Modulation now was measured at dis- 
tances of 300, 900, 1500, and 2100 feet by moving the receiver laterally 
in 18-inch increments, placing it in the four appropriate positions to 
intercept the beam. 
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Fig. 6— Range dependence: multiple pass (300-2100 ft); 
(a) schematic, 
(b) data 8/6/65 night—8 runs, 
8/8/65 night—8 runs, 
8/12/65 night—5 runs, 
8/17/65 night—10 runs. 
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On some nights, beam dancing was severe enough to carry the beam 
off the final mirrors periodically, so that data could not then be taken. 
In the early part of most summer nights there was a slow drift of 
the beam upward as the air cooled increasingly below the temperature 
of the roof and the strength of the inverted “prism” increased with 
time (the situation is different during winter nights, see Fig. 11). A 
representative speed of vertical beam movement at 2100 feet was about 
two inches per hour. This required us not only to align the mirror 
system for each night’s work but to realign each night every 20 to 30 
minutes. Mechanical stability of the mirrors was good. When atmos- 
pheric conditions were stable enough there was no noticeable move- 
ment of the beam, so that building vibrations appeared to have no 
important effect. 

The averages of the data were plotted in Fig. 6(b). It is seen at 
once that there is no significant change in slope beyond 300 feet. The 
points at 300, 900, and 1500 feet yield a well-defined straight line 
of slope 1.39, in good agreement with results below 300 feet. The value 
at 2100 feet seemed distinctly low, and it was not used in fitting the 
data. This is reminiscent of the behavior of the last point in Fig. 3 
which led to the present measurements. It was not possible to deter- 
mine any systematic cause of error in the measurements at the 
largest distance. The fact that there the beam is largest, and most 
difficult. to collect, would explain a modulation which is too large, not 
too small. 

Once again an increase in the distance seemed necessary, and this 
was accomplished by adding one more 300-foot pass to the folded 
system and moving the receiver to the transmitter end. This is shown 
in Fig. 7(a). Here the first three mirrors were four inches square, and 
the last four were six inches in diameter. The data, averaged in the 
usual way, are presented in Fig. 7(b). The slope of 1.46 again agrees 
with previous values, within experimental error. The value at 2400 
feet shows that there is no decrease in slope beyond 2100 feet. It 
would seem that the 3/2 power dependence of modulation depth on 
distance continues to at least 2400 feet. It would be instructive to 
continue these measurements at larger distances, but, for this, larger 
optics would be required than were available. 

The folded-path method of obtaining large ranges is both convenient 
when space is limited and desirable when readings must be made 
quickly at positions which would be widely spaced in a straight path. 
However, the question remains whether this is in every sense equiva- 
lent to the unfolded, straight path. For the distances considered here, 
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Fig. 7-—- Range dependence: multiple pass (600-2400 ft); 
(a) schematic, 
(b) data 8/18/65 night—8 runs, 
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it takes only about 0.3 microsecond for the light to travel each lap 
and thus a total of 2.4 microseconds for 8 laps of 2400 feet. The char- 
acteristic time of turbulence is of the order of milliseconds, and con- 
sequently the turbulence can be considered frozen during the time 
required for the light to travel the 2400-foot folded path. Under the 
circumstances, not only are there possible correlations between the 
closely adjoining segments of the folded path but, since the beam 
travels oppositely along the successive segments, there may be an 
actual reduction in the net effect of atmospheric turbulence. 

To determine this, we now set up a mutiply-folded system of mir- 
rors, in all respects the same as the 2100-foot system in Fig. 6(a) 
except that now the segment length was reduced from 300 feet to 38 
feet. This gave a total, folded distance of 270 feet so that comparison 
could be made with the results for a straight 300-foot path which 
were given in Fig. 3. The new arrangement is shown in Fig. 8(a). The 
results appear in Fig. 8(b). The functional dependence is sensibly 
unchanged from that in Fig. 3. This therefore not only vindicates the 
folded-path method used for large distances, but it further confirms 
the short-distance behavior of the modulation. 

Although this experiment provided no evidence of correlation effects, 
it seemed worthwhile to make a further attempt to detect any correla- 
tion effects that may affect depth of modulation. The folded path, 
therefore, was rearranged as shown in Fig. 9(a). The telescope was 
focused on the receiver to give the smallest possible spot, which was 
about one cm in diameter. The first, second, fifth and sixth mirrors 
were two inches in diameter and quarter-wave flat. They were 
mounted so that no frame of the mirror or other obstruction extended 
in front of the six-inch diameter third and fourth mirrors. It was thus 
possible to position the spots on each bank of mirrors just two inches 
apart between centers. Hence, the average separation of successive seg- 
ments of the path was only one inch. However, this made it impossible 
to intercept the beam with the photomultiplier, as was done before, 
without obstructing a previous segment of the beam. Consequently, 
the beam was reflected to the photomultiplier with an elliptical mirror 
(of the telescope-diagonal type). This mirror was then displaced 
laterally to intercept the beam. In this experiment, the positioning of 
the mirror was simple enough that a run could be completed in about 
2 minutes, so that many runs were quickly made and exceptionally 
good averaging should result. 

The average values are plotted in Fig. 9(b). Once again, the now 
familiar features of the curve appear. It is interesting to note that, in 
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Fig. 8— Range dependence: multiple pass (42-270 ft); 
(a) schematic, 
(b) data 10/7/65 day—d runs, 
10/13/65 night—10 runs, 
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the previous experiments, the slope beyond 100 feet has centered 
around 3/2 and that, in the present experiment in which especially 
good averaging against changing conditions was expected, the slope 
was 3/2 as closely as such curve-drawing will allow. The present re- 
sults again give no evidence of correlation effects. 
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Fig. 9 Range dependence: multiple pass (42-270 ft) with mirror separation of 2’; 


(a) schematic, 

(b) data 8/27/65 night—7 runs, 
8/30/65 night—11 runs, 
9/2/65 night—13 runs. 
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V. RESULTS—DEPENDENCE ON ATMOSPHERIC CONDITIONS 


It was remarked before that the first and dominating result of any 
measurement of depth of modulation is that the value changes steadily 
with changing atmospheric conditions. This persists during, and inter- 
feres with, all attempts to measure the dependence of modulation on 
any other parameter. Hence, a great deal of qualitative observation of 
the effects of conditions was inescapably made during the measure- 
ments of aperture and range dependence. And this extensive experi- 
ence taught simply that the depth of modulation was unpredictable. 
There never appeared any simple correlation between modulation and 
qualitative observables such as wind speed or direction, sun conditions, 
or time of day. The most promising lead which developed was the 
observation already cited in the aperature work that the depth of 
modulation was markedly smaller over a 25-foot range after a short 
heavy rain than before [Fig. 2(b)]. The rain, of course, would have 
both cooled the roof (which had been heated by the sun) and reduced 
dust in the air. Since no connection previously was apparent between 
modulation and temperature conditions, it then was thought that the 
modulation might be affected strongly by particulate scattering. 

Indeed, at night the beam was always well decorated by forward 
scattering from what appeared to be dust and haze particles. Normally 
the beam could be seen with the eye placed within about 10-* radian 
of the forward direction. This was the case with clear or hazy con- 
ditions. Only with fog could the beam be seen substantially more than 
10-? radian away from the forward direction. Under no conditions, 
including light fog, could backscattering be detected by eye. The 
result with the multi-folded system was that, when standing at one 
end of the range and looking toward the other, only those segments 
in which the beam was approaching could be seen. This is shown by 
the photograph of Fig. 10, which was taken looking toward the trans- 
mitter. Speed of the film was ASA 125. Exposure time was 15 minutes, 
so the beams appear well filled-in. There seem to be four separate, 
parallel beams from four separate sources. There is no trace of the 
diagonal connections of the returning beams. It seems clear that back 
scattering was orders of magnitude weaker than forward scattering. 
The right-hand beam came from the laser, and light scattered at the 
source caused the saturation of the film. 

When measuring depth of modulation, therefore, the level of for- 
ward scattering out of the beam was readily observed. Attempts were 
made to correlate modulation with scattering, and no qualitatively 
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Fig. 10— Photograph of multiple-pass laser beam illustrating the lack of back 
scattering. 


apparent relation could be found. The level of scattering did vary in 
a pronounced way as the density of particles varied, but modulation 
did not seem to change correspondingly. With short ranges, such as 
25 feet, cigarette smoke was deliberately blown into and along several 
feet of the beam causing brief and strong decoration of the beam but 
having no noticeable effect on modulation. 

Having thus found no dependence by casual observation, we ar- 
ranged to measure modulation while also systematically noting wind 
speed, direction, and variability, temperature of the roof, and temper- 
ature of the air at beam level. The arrangement was that shown in 
Fig. 3 except the distance was 138 feet. Temperatures were measured 
at the middle of path. Readings were taken every half hour continu- 
ously for 24 hours so that any diurnal variation would be detected. In 
particular, one might expect a change in modulation depth around sun- 
rise and sunset. (Reliable observations of “good-seeing” at sunrise and 
sunset date at least from 1878 when Michelson measured the velocity 
of light and found that he could not work at any other time due to 
excessive “boiling” of the image of his slit.§) 

The 24-hour period began at 5:45 p.m. on November 3, 1965. A 
broad variety of wind conditions was obtained both at night and dur- 
ing the day, ranging from dead calm to a period which was violently 
gusty in the late morning. The sky was clear at the beginning, becoming 
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gradually overcast until no stars could be seen by 3:45 a.m. on No- 
vember 4, 1965. This cleared to only a hazy sky by 8:00 a.m., and the 
daytime then remained clear except for occasional clouds. 

The results are assembled in Fig. 11. The difference between roof 
and air temperature is recorded since one expects temperature gra- 
dients to be more significant than temperature level. Wind conditions 
are given below depth of modulation. It is apparent at once that there 
was no simple correlation between modulation and temperature diif- 
ference or time of day. Any distinct reduction or other effect on mod- 
ulation at sunrise or sunset is conspicuously absent. 

The negative results of the previously qualitative observations 
therefore have been extended by these more quantitative results. It is 
still clear that modulation depth depends upon atmospheric conditions 
in a sensitive way, but the measurements so far have not revealed the 
nature of the dependence. 


VI. THEORETICAL BACKGROUND: DISTANCE DEPENDENCE 


The amplitude and phase fluctuations of an electromagnetic signal 
that has propagated through a random medium have been theoretically 
analyzed using different approaches. Of the more familiar ones are the 
following: the ray theory, the first Born approximation, and Rytov’s 
method. The results obtained by all three methods agree, though they 
may hold true only for certain regions of distance.® Of the above 
three approaches, Rytov’s method lends itself to a generalized treat- 
ment that holds good for both short (Fresnel zone) and long (Fraun- 
hofer zone) distance regions of the scatterer. It should be noted that in 
the following discussion the Fresne] and Fraunhofer regions are with 
respect to the scattering medium and not with respect to the trans- 
mitter aperature. 

The physical configuration for which the calculations are made is 
shown in Fig. 12. An infinite plane wave is incident upon a semi-infinite 
(-~x7 <4 < +03; —-—a~ < y < +0;2 = 0) inhomogeneous and ran- 
dom medium which is assumed to be quasistatic. A point detector 
is located at « = L which not only sees the unscattered direct wave 
but also waves scattered from within a scattering volume that is in 
the form of a cone. This cone has its vertex at the receiver and has 
an aperture angle of the order of 1/ka, where k is the propagation 
constant (= 27/)) and a is the scale size of the tubulence. The justi- 
fication of using this configuration for our measurements will be given 
at the end of this section. 
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Fig. 12— Physical configuration for calculation of amplitude fluctuations. 


The functional dependence of the amplitude fluctuations will depend 
on whether the receiver is located in the Fresnel or Fraunhofer zone. 
The extent of these zones are determined by a dimensionless parameter 
D, called the wave parameter. It is given by 


4D, 
~ ka? _ 
D <« 1 for the Fresnel region and >1 for the Fraunhofer region. The 
dependence of the mean square fluctuation of the amplitude B’ on this 
distance LZ is given by” 


For D<1 Bal 
D>1 Bal. 


In other words, the root mean square value of the fluctuations will 
have a functional dependence on distance given by the distance raised 
to the power 3/2 and 1/2 for the Fresnel and Fraunhofer zones, re- 
spectively. 

The quantity measured in the experiment is the percentage modula- 
tion given by 100 times the ratio of the rms voltage to the average 
value. This ratio is the same as the ratio of the rms value of the fluctu- 
ating light power to the average light power expressed in percentage. 
For values of ac-to-de power ratio very small, 1t can be shown that 


Vous _ ere _ fe oie 
= eed 


Ves i ee Lay (3) 
where V’s refer to the voltages across the photomultiplier load, P’s to 
the light power and £’s the electric field intensities in the light radia- 
tion field. From (2) and (3), we see that the ratio of the Vims to Vay 
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should vary as L°/? in the Fresnel zone and as L!/” in the Fraunhofer 
zone. Even though the configuration used for the theoretical calculation 
does not seem to represent our experimental arrangement, its validity 
can be justified by the following argument. The 3/2 power law, which 
is the one that will be used to compare with our results, can also be de- 
rived using ray theory in which fluctuations in an infinitely thin ray of 
light are considered. The fluctuations in this case can be measured us- 
ing a point detector. In our case, the ray is of finite diameter (that 
corresponding to the laser beam diameter), and consequently the detec- 
tor is also of finite dimension to insure collection of the entire ray. 


VII. DISCUSSION 


The modulation of a single-mode, single-frequency laser beam at 
6328A° by atmospheric turbulence was investigated by varying the 
propagation distance as the parameter. Unlike the spectral width of 
modulation, the depth of modulation does depend on distance and 
varies as 3/2 power of the distance. In making the measurement, it 
was ensured that all of the direct beam was collected by the receiver. 
The range of the propagation distance was extended from 0 to 2400 
feet. The empirically obtained 3/2 power law agreces well with the 
theoretical result obtained using Rytov’s method, provided the propa- 
gation distance is within the Fresnel zone of the scatterer. This assump- 
tion leads to an estimation of the scale size of the atmosphcric turbu- 
lence. The effective scale size in (1) is estimated to be larger than 5 
cm in diameter for L = 2400 feet and D = 0.1. 

An interesting feature of the dependence of depth of modulation on 
distance is its large value at distances lower than about 100 feet. This 
produces a functional dependence on distance which is other than 3/2 
power. We have also noticed that the spectral width which is char- 
acteristicaly a few hundred hertz and whose exponential decay with 
frequency is otherwise independent of distance undergo changes at 
these short distances. The amplitudes of the low-frequency compo- 
nents decrease more rapidly than those of the high-frequency com- 
ponents as the distance is made shorter. We believe that the short 
distance variation of the depth of modulation and the spectral width 
are related and are caused possibly by the same phenomenon. This 
is under further study. 

In comparing the absolute value of depth of modulation with those 
obtained by others, it has to be borne in mind that we collect all of 
the direct beam in contrast to previous work in which only part of the 
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direct beam was collected. Consequently, the magnitude of the depth 
of modulation measured by us is much less than that of others. For 
example, Edwards and Steen! have observed with a zirconium arc 
source for a 300-meter path a depth of modulation as high as 80 
percent, whereas our value for the same distance is of the order of 
1 percent or less. Another measurement by Portman, et al’® with a 
partially collected beam at 500 meters yielded a peak to peak percent- 
age modulation of 150 percent. 

Although the main goal of this work was the functional dependence 
of the depth of modulation on distance, two other observations were 
made which are worth noting. Contrary to the behavior of the spec- 
tral width, which depends systematically on weather conditions, the 
depth of modulation has no clear-cut dependence on the atmospheric 
conditions which were measured so far. It seems modulation is a sensi- 
tive but obscure function of atmospheric conditions. The atmospheric 
variables which need to be measured evidently are relatively fine. 
Besides many qualitative observations, the quantitative results of a 
24-hour run were a demonstration of this. 

Also, we have not been able to observe any back scattering of the 
laser radiation even under severe weather conditions. Most of our 
experiments were conducted during night time. Even on very dark 
nights, the dark adapted eye (of several observers) could not detect 
any trace of back scattering. This is true in conditions of clear atmos- 
phere with various amounts of particulate matter, haze, fog, and 
under severe rain storms. From these qualitative observations, we 
are led to estimate that the back scattering is orders of magnitude 
lower than the narrow angle forward scattering—a value considerably 
lower than the 2 percent obtained by Carrier and Nugent.!* This 
observation is surprising also in view of the various reports of at- 
mospheric back scattering observed with optical radar systems (e.g., 
Collis and Ligda?*). 
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Iterative Solution of Waveguide 
Discontinuity Problems 


By W. J. COLE, E. R. NAGELBERG and C. M. NAGEL 
(Manuscript received November 25, 1966) 


The application of matrix iterative analysis to the solution of waveguide 
discontinuity problems is discussed. It 1s concluded that the ‘Gauss-Seidel’ 
or “‘noint-single-step’’ method offers several advantages over more conven- 
tional invertive procedures, particularly in the speed of execution. Two 
examples are presented as illustrations: analysis of an H-plane discon- 
tinuity an a rectangular waveguide and conversion from TE,, to TM,, 
modes at an abrupt discontinuity in a circular waveguide. The latter 
results are shown to be in good agreement with measured values obtained 
in a previous investigation. 


Y. INTRODUCTION 


The analysis of waveguide discontinuities, for application to the 
design of antennas and microwave networks, continues to offer challeng- 
ing problems in electromagnetic theory and microwave engineering. 
Thus far, the solution of these problems has depended to a large extent 
on various approximate techniques, such as variational and quasi- 
static methods,’ which are extremely useful but nevertheless limited 
in applicability. 

The shortcomings of classical analysis have been surmounted to a 
large extent by our ability to solve electromagnetic boundary value 
problems by numerical methods, making extensive use of digital com- 
puters. Computational techniques are not only an abundant source 
of engineering data, which might otherwise require elaborate construc- 
tion and experiment, but they can also provide a unique analytical 
laboratory in which to evaluate approximate theoretical methods under 
easily controlled conditions. In this paper, we shall be concerned with 
these numerical methods as they apply to certain waveguide discon- 
tinuity problems. 
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For the sake of simplicity we shall consider, as an example, the 
problem of two waveguides with similar cross sections connected to- 
gether at the plane z = Q, as illustrated in Fig. 1. A wave is shown 
incident from the smaller waveguide impinging on the discontinuity. 
The result will, of course, be to excite an infinite number of normal 
modes in each guide, some of which carry real power away from the 
junction, with the remainder being evanescent and contributing to 
the electromagnetic field only in the vicinity of the connecting aperture. 
It must be recognized that these evanescent modes play an important 
role since they, in part, determine the amplitudes and phases of the 
propagating modes. It is the fact that an infinite number of waves 
must, in principle, be considered that makes this type of problem so 
difficult. 

The contents of the paper may be summarized as follows: We begin 
by establishing an appropriate form of the uniqueness theorem for 
Maxwell’s equations as they apply to boundary value problems of 
this type. In numerical analysis, the criteria for uniqueness are of more 
than academic interest since they provide meaningful and practical 
methods by which to assess the accuracy of results. Next, the normal 
mode representation of the fields is discussed, the object being to 
arrive at a matrix equation formulation of the problem in which the 
components of the unknown vector are the modal coefficients. It is 
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Fig. 1—- Waveguides of similar cross section connected at the plane z = 0 by 
an abrupt discontinuity. A wave is assumed incident from the smaller guide. 
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suggested that this matrix equation may be solved by an iterative 
procedure and, upon studying the convergence properties of such 
methods we find a critical dependence on the particular algorithm 
used. Two examples will be presented as illustrations, analysis of an 
H-plane discontinuity in a rectangular waveguide, and conversion from 
TE,, to TM,, modes at an abrupt discontinuity in a circular waveguide. 
The latter results are shown to be in good agreement with measured 
values obtained in a previous investigation. 

Rationalized MKS units and the (suppressed) harmonic time de- 
pendence exp (—2wt) will be used, unless otherwise specified. 


II. UNIQUENESS AND ERROR CRITERIA 


A representation of the discontinuity is shown in Fig. 2. It is assumed 
that the regions to the left (denoted by —) and to the right (denoted 
by ++) are each filled with homogeneous material, but with possibly 
different constitutive parameters. Maxwell’s curl equations in the 
respective regions are thus given by 





(1) 
+ . 0E* 
V XH’ =e ey 


As usual for uniqueness theorems, we begin with two solutions in 
each region presumed to be correct, and denote the differences respec- 
tively by E*, H*.* Then from the Poynting theorem,” it follows that 





Fig. 2— Waveguide discontinuity showing boundary surfaces A, B, C, D and 
respective normals ni, Me, Ms, M4. 


_* Physically, these fields would correspond to a waveguide discontinuity problem 
without excitation. 
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in the — region 


If (E, X H,):n, dS + If (E; X H>)-n dS 


_ akin tes 2eHd apiece 
ee vill. 2-2 Wok va dv (2) 


and in the + region that 


If (EY X Hi)-n, dS + If (Et X H*)-n, dS 


£0 + + Pots. + 0 + yt 
~ valle: ee aes oe 


where 0/dt denotes differentiation with respect to time and the sub- 
script ¢ denotes the field transverse to the generatrix of the cylinder. 
The unit normal vectors n, , n., n; , and n, are shown in Fig. 2. 

One must also take into account the fact that certain physical con- 
siderations will limit the class of admissible solutions. For example, 
if we let the surfaces A and D recede to infinity, then all evanescent 
modes will have decayed to zero and the respective surface integrals 
then represent the power flow away from the discontinuity. Assuming 
no loss, the total power must vanish. Furthermore, it can be shown 
from Maxwell’s equations that the transverse components of electric 
and magnetic field at the interface must be continuous. Adding (2) 
and (3), we find that the following time derivative must vanish, 


ale IH. EE dV+yue IH. H’-H™ dV 
ig If. En dv ey? [[[ wax av | -0, (4) 


We may, however, regard the quantity in brackets as having had 
a zero value at some time, say at ¢ = 0, the excitation time. The term 
in brackets therefore, vanishes for all time and, since each of the in- 
tegrands is positive semi-definite, they must vanish separately. Thus, 
at each point in the + and — regions, 


E; — E; = Hi — H, =0 
E, — E,; = H,; — H, = 0 


(5) 


and the solution is thereby shown to be unique. We may now state 
the following uniqueness theorem for waveguide discontinuity problems. 
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Theorem: The solution to a waveguide discontinuity problem is uniquely 
specified if 7¢ can be shown to have the following properties: 


(2) Lt satisfies Maxwell’s equations and the appropriate boundary 
conditions in the regions on each side of the discontinuity. 
(1t) The components of electric and magnetic fields tangent to the 
interface are continuous. 7 
(zit) In the case of a lossless discontinuity, energy 1s conserved. 


These three conditions obviously play an important theoretical role 
in the solution; where numerical methods are used, they also provide 
fundamental criteria by which the accuracy of computed results can 
be assessed. Accordingly, we shall define the following quantities to 
be used as error criteria: First, there is the parameter €p , which indicates 
how well the solution conserves energy, given by 


Pent. 
PP ins 


where P, , P,, and P;,, are the reflected, transmitted, and incident powers, 
respectively. Second, the mean square error in the tangential electric 
field is defined by 


I, (6) 


Cp = 


[[  |\@-By Pea 
ey = Aperture 
I[ | EB i dA 
Aperture 


and third, for the magnetic field, 


|| | Hin”? P dA 
t 
Aperture 


where E“"” and H"“"” refer to the incident wave. 

The smaller the quantities ep, €,, and e€,, the more closely the 
boundary conditions are satisfied, at least in the mean square sense, 
and the more accurate we shall consider the solution to be. 


(7) 


Ill. MATRIX FORMULATION OF THE BOUNDARY VALUE PROBLEM 


The most convenient format for numerical solution of waveguide 
discontinuity problems is a matrix representation, in which the modal 
coefficients form the unknown column vectors and the discontinuity 
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is characterized by a square matrix. We recognize that there is also 
an analogous integral equation in terms of the aperture electric or 
magnetic field. However, since the numerical solution of the integral 
equation is generally carried out by reducing it to a matrix equation, 
we shall proceed to the matrix formulation directly from the physical 
characterization of the boundary value problem. This matrix equation 
will then be solved by an iterative method, the theory of which is 
discussed in Section IV. 

It is assumed that in each of the waveguides, the electromagnetic 
fields may be characterized by a denumerable set of known vector 
eigenfunctions which may be ordered according to some index. We shall 
be concerned only with the transverse fields,* denoted as follows: 


"E’(r) (p = 1, 2, 3, ---) denotes the transverse electric field for the 
pth TM mode in the + waveguide, with r as the position vector in 
the transverse plane. 

"H/(r) = transverse magnetic field for the pth TM mode in the 
-+ waveguide. 

*E/’(r) = transverse electric field for the pth TE mode in the + 
waveguide. 

“H/’/(r) = transverse magnetic field for the pth TIX mode in the 
+ waveguide. 


By replacing the + by — we have the analogous notation for the 
other waveguide. An important point concerning sign convention is 
that the unknown modes in the — waveguide will all be taken to 
propagate away from the discontinuity, l.e., in the —z direction. 
Although the electric field does not change sign when the direction 
of propagation is reversed, the magnetic field does, and this fact must 
be carefully taken into account. 

In order to define the amplitudes of the respective vector wave 
functions, we adopt the following normalization,” written in terms of 
integrals over the waveguide cross sections: 


i “EB dA = | *h PF dy, (9) 
tA 


i “Ey! *Ei/* dA = wp’ bya (10) 
=A 


in which A, is the respective characteristic wavenumber, and uy is the 
permeability, which in our case will be the permeability of vacuum, 


* Tt is assumed that the individual waveguides can support pure TE and TM 
modes, which is the case for applications of interest here. 
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since both waveguides will be assumed empty.f By introducing the 
Kronecker delta 6,,, we have also expressed the fact that the trans- 
verse fields in the individual waveguides are orthogonal. 

Once the normalizations for the electric wave functions are defined, 
those for the magnetic field are also specified since, for both the TI 
and TM modes, the transverse electric and magnetic fields are uniquely 
related. In particular, for a TM mode 





*H’ = —S e, X *E! (11) 
and for a TE mode 
‘HY = + hy e, X “Es’, (12) 
COU 


where e, is a unit vector in the z direction. Note again that the sign 
convention 1s such that a field in the — waveguide is taken to be a 
reflected wave, travelling away from the discontinuity. The magnetic 
field normalization is thus given by 


i “W.*H dA = we 5, (13) 
=A 


/ “HY. tH dA = | tA Pb, (14) 
+A 


Both sets of transverse wave functions have the property of com- 
pleteness, which is to say that any transverse electric (or magnetic) 
field can be synthesized from a set of TE and TM vector wave functions, 
provided that the directions of propagation of the normal modes are 
known. For the problems to be considered here, this latter information 
is available from physical considerations, since all modes propagate 
away from the Junction with the exception of the incident wave whose 
amplitude is known. This amplitude will be taken to be that of a 
normalized mode. 

We now derive the appropriate matrix representation for the dis- 
continuity problem. Assume a dominant (TIi) mode wave (E/’, H/’) 
is incident from the — guide, setting up a transverse electric field in 
the aperture just to the left of the junction. This field, referred to as “E, , 
may be synthesized as follows: 


“E, = “E’ + 3074, EJ + 3S BSE!’ (15a) 
pal qg=1 


+ The asterisk (*) denotes the complex conjugate. 
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with the modal coefficients ~A, and ~B, as yet undetermined. The 
corresponding transverse electric field on the + side, denoted by “E 
would then be given in terms of normal modes on the + side by 


*E, = 014, E, + DS BE’. (15b) 
p=1 q=1 
Since the transverse electric field is continuous across the aperture, 
we have that 


"E, = E, on C 


"E, =0 on D-C. 


As shown in Fig. 2, D-C represents the conducting wall which makes 
up the remainder of the junction, and on which the transverse electric 
field must vanish. Expanding (15a) in a Fourier series of modes in 
the + waveguide, we find that the modal coefficients are related by 


"A, = TRE we |. “E/’. (ES dA + Tara ; i a a “E/. "El* dA 





Bite =» B, [ ECE dA (16a) 
7 C 


Dp 


fie i tao : “ER dA a A, | “BEI dA 
@ Uh C Wh gq=1 Cc 





+= 7B, | “EY “E*dA — (16b) 
C 


Ww pL q=1 


or, More succinctly, using partitioned matrix representations, 


*® 1 Fo 











0 —z-35 J 
OM 
ef ™y 4 ae ' ’ + 
+h! wa / / vr / = 
+ ihe (17) 
7 5 J — g et ’ + er ae + B 
@) ph row noo 
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in which the vectors and submatrices are defined as follows: 
As *B, 
“@ ome "Aas *B — *B (18) 


(5,), = | “BYE dA 
Cc 





(19) 
(Gin. | “BYE! dA 
Cc 
1 
Pap OS” 
+ , 1 
D = 0 | thi |? 0 ma (20) 


and g is the identity matrix. The matrix &, whose transpose appears 
in (17), is the matrix of coupling coefficients, defined as the scalar 
products of electric transverse vector wave functions for the wave- 
guides on each side of the discontinuity. The four index notation Js 
interpreted as: 


; 4? 
? 


Sng 7 7 = | “EB! . *E//* dA (21) 
Cc 


with analogous definitions for other combinations. 

The system of equations given in (17) is clearly underdetermined 
since the number of unknowns is twice the number of equations. How- 
ever, an additional set can be derived by employing the boundary 
condition that the transverse magnetic field must also be continuous 
across the interface. The matrix equation, analogous to (17), but 
corresponding to this second boundary condition, is given by 


1 ] +,- ae(te 

-@ 0 oe ge ne vt J @t 

| a : : (22) 
Ol LG) Pa op! eft eh] @* 


4 





4} 4} 
? 
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in which 
1 
0 
ies (23) 
0 
1 
i nn 
Saye) Gi tees: 0 doe (24) 
| ha" | 
0 0 


and the matrix 3C, whose transpose appears in (22), is the matrix of 
scalar products of magnetic transverse vector wave functions. The 
four-index notation is interpreted in the same way as in (21). 

It should be noted once again that all the matrices which appear 
in (17) and (22) are infinite matrices, corresponding to the fact that 
in general an infinite number of modes are excited in the neighborhood 
of the discontinuity. In practice, of course, there must be a truncation 
and the problem then becomes one of solving a set of matrix equations 
whose order, NV, depends on the accuracy required. Unfortunately 
there is, as yet, no way in which the number of modes required to produce 
a given accuracy can be predicted. We can only emphasize the need 
for meaningful error criteria which will act as a guide in choosing a 
number of modes which will be large enough to give sufficiently accurate 
results but at the same time not be so large as to require excess com- 
putation. It is expected that the criteria given in Section IJ will prove 
very useful in this respect. 


IV. MATRIX ITERATIVE METHODS 


It was shown in the previous section that the waveguide discontinuity 
problem of interest here can be formulated in terms of a system of 
linear algebraic equations. This is a recurrent theme in mathematical 
physics, so that an extensive theory concerned with the efficient solution 
of matrix equations has evolved. In this section, we shall be concerned 
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with some of the elements of this theory, placing particular emphasis 
on the solution of matrix equations by iteration. 

The system of linear algebraic equations which results from satisfying 
the aperture boundary conditions on the transverse electric and mag- 
netic fields can be written in the matrix forms 


“a =U+t+R-a (25) 
“a =U+ 8-"a, (26) 
where 
| *@ 


= | | si 


the vectors U and VU and the matrices ® and § being correspondingly 
identified from (17) and (22). Equations (25) and (26) are easily un- 
coupled to give 


[J — RS]-"a = U+ KV (28) 
[Ig — SQ]--a@ =OU+ Su (29) 

both of which are seen to have the general form 
Mae = y. (30) 


In (80), x is an N-dimensional complex vector whose components 
are the coefficients of the normal modes in the two waveguides, 9W is 
an N X N complex matrix characterizing the discontinuity, and y is an 
excitation vector duc to the incident wave. 

The obvious method of solving (80) is to compute the inverse of 3% 
and thus directly obtain 


x= My. (31) 


However, we should recognize that it is the solution vector x which 
is required, and that computing the inverse is not always the best equa- 
tion solving technique. For example, because the modal coefficients 
may decrease slowly with mode index, an accurate approximation of 
the physical problem often requires that 9% be a very large matrix, 
and inversion procedures for large complex matrices require considerable 
computational effort. An alternate approach is therefore suggested, 
namely the solution of (80) by a method of iteration. 

In an iterative algorithm, we begin with an initial ‘“‘guess” for the 
solution and, from this, generate a supposedly improved solution, 
repeating the process until successive iterations give results which 
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we. 
rae 


~*TE,9 MODE INCIDENT 


Fig. 3—H-plane discontinuity in a rectangular waveguide. The incident mode 
is &@ TE) mode (electric field vertical). 


agree to within some prescribed norm. The solution to which the 
procedure converges must, of course, be independent of the initial as- 
sumption. 

A tempting iterative procedure for the present problem is suggested 
by writing (28) in the form 


‘Ta = U+ RV 4+ RS8:*a (32) 
with an initial assumption 
‘a = Wt RV. (33) 


Physically this corresponds to first assuming the aperture electric 
field to be that of the unperturbed incident wave, and calculating the 


REAL PART OF *TE\9 COEFFICIENT 





NUMBER OF ITERATIONS 


Fig. 4— Transmission coefficient of an H-plane discontinuity. 
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TE and TM modal coefficients in the + waveguide on this basis. The 
corresponding magnetic field is then determined on the + side of the 
aperture and, with the aid of the appropriate continuity condition, 
is used to find the magnetic field and subsequently an “improved” 
electric field in the — region. This second guess for the aperture 
electric field is then used to repeat the process, etc. 

As a test, this algorithm was applied to analysis of an H-plane dis- 
continuity in a rectangular waveguide, the incident wave being a 
TE,. mode of normalized amplitude [see (14)] as in Fig. 3. The dimen- 
sions were ka, = 4.5 and ka, = 3.5 where k is the free space wave- 
number. With this choice of parameters, only the TE,) modes can 
propagate in each guide. (This problem is discussed in further detail 
in Section V.) 

Fig. 4 shows the result of calculating the real part of the modal 
coefficient for the TE,, mode in the larger waveguide, plotted as a 
distribution of points giving the value at each iteration. The Fourier 
series for this particular example was truncated after twenty five terms. 
Aside from a small amplitude oscillation of less than one percent rms, 
the results seem reasonable, especially in view of calculations for the 
mean-square errors €y and ¢,, which are illustrated in Fig. 5. These 


0.032 


0.012 


MEAN SQUARE ERROR 





NUMBER OF ITERATIONS 


Fig. 5— Mean square errors in transverse electric and magnetic fields for an 
H-plane discontinuity. 
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REAL PART OF *TE,g COEFFICIENT 





NUMBER OF ITERATIONS 


Fig. 6— Oscillatory instability of a nonconvergent algorithm for the H-plane 
discontinuity analysis. 


decrease monotonically with succeeding iterations, ¢, approaching 
approximately 0.001 and ey, a value of about 0.012. The larger error 
in the magnetic field can be attributed to the singular behavior in 
H near the corner of the discontinuity. The asymptotic value of the 
energy parameter €p is approximately 0.007. 

The apparently accurate results obtained using this algorithm are 
in fact quite deceptive, and may actually be attributed to the propitious 
initial choice for the aperture electric field. It will be recalled that a 
very important criterion for validity of an iterative procedure is that the 
results be independent of the initial assumption. In order to determine 
whether such a criterion is satisfied for this particular algorithm, the 
TE. modal coefficient for the larger waveguide was arbitrarily doubled 
after the fourth iteration, which is equivalent to deliberately assuming a 
poor initial choice for the aperture field. The effect, shown in Fig. 6, 
indicates that the algorithm does not relax to the previous values, but 
continues to oscillate with a large amplitude. Similarly large fluctua- 
tions occur in €g , €,, and ép, the conclusion being that this particular 
precedure is not satisfactory, and can be expected to give reasonable 
results only if the initial choice is a very good one. The reason for this 
instability will become apparent after we consider those aspects of 
matrix-iterative analysis which are relevant to these problems. 

In the usual framework for iterative procedures, the matrix equation 
satisfied by the unknown vector x is written in the form 


where 3 and f are appropriate to the particular scheme being used. 
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This leads very naturally to the recursive formula relating the m + 1 
to the mm iteration, 


gt) = ote + f. (35) 
We denote the error vector at any iteration by e‘”’, where 
e™ = o™ — Z, (36) 


and € is the exact solution to (84). Then, by substituting (86) into (35), 
we find that e‘”*” is related to e“” by 


eo Mite. (37) 


Therefore, the error at the mth iteration is expressible in terms of the 
initial error by 


er SS oe. : (37a) 


For an absolutely converging solution we thus require that 


e™ 30 as moo (38) 


(0) 


regardless of the initial guess x°’. This is equivalent to 


wm” —-O0 a mao (39) 


where the 0 in (88) and (89) denotes a null vector or matrix, respectively. 
It can be shown* that an N X N complex matrix 9M is ‘convergent’, in 
the sense of (39), if and only if all the eigenvalues A, of 1 magnitude 
less than unity, i.e., 


ele at: 2. (40) 


We can easily see why this requirement will guarantee convergence, 
at least for the special case where the eigenvectors a; of SW span the 
space of N-dimensional complex vectors. The initial error is then 
expressible as 


N 
a D Cia; , (41) 
i=l 


where the C’; are constants. The error at the mth iteration then be- 
comes, from (37), 


N N 
am = >» CM" a, = > C Nia; (42) 
i=l i=1 


and, as m —> o, each term in the sum approaches zero, provided, 
of course, that (40) is satisfied. 
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Equation (42) also reveals useful information concerning the speed 
of convergence, which is seen to depend on how close the magnitudes 
of the A; are to unity. Clearly if the largest eigenvalue is very near 
one in magnitude, the convergence will be very slow and hence a large 
number of iterations will be required. It would be logical, in the light 
of this reasoning, to evaluate the magnitude of the largest eigenvalue 
for the H-plane discontinuity problem discussed proviously. Unfor- 
tunately, because of the geometrical asymmetry, the matrix 9% 1s not 
Hermitian and so the usual computational techniques for determining 
eigenvalues cannot be used. We can, however, find an upper bound 
for the modulus of the maximum eigenvalue, p(9%) = max {| , |}, 
given by” 


(Mt) S [e(an'sn)]}’. (43) 


where { denotes the conjugate transpose matrix. Note that MN'9 is 
Hermitian, so that standard computer programs can be used to evaluate 
its eigenvalues. We find for the previous H-plane problem that | A; | S 1.16 
which, although, not conclusive, shows the possibility of such an os- 
cillatory instability. 

One technique which is suggested as a means of obtaining a con- 
vergent algorithm is called the “Richardson” or “‘point-Jacobi” method.” 
In this approach, the matrix 9M is first partitioned as 


M=O+L4+ 4, (44) 


where D is a diagonal matrix containing the diagonal terms of I, £ 
is a strictly lower triangular matrix and U is a strictly upper triangular 
matrix. The system (34) is then written as 


(§- De=(L+ We ti (45) 
from which 
(§— D) (Lt We + (§ — DF 
Wer + fr, (46) 


& 
| 


I 


(46) being the matrix representation of the “Richardson” or “point- 
Jacobi” method. 

A modification of this procedure is referred to as the “Gauss-Seidel” 
or “‘point-single-step” iteration method.° Note that in solving (46) by 
iteration, the components of 7‘"*” are all computed from the com- 
ponents of x‘. Intuitively, it would seem more attractive to use the 


latest estimates of z, 1e., In computing 2°”*’ we should use, wherever 
P b] 7 ? 
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they appear, the components 2;"*" (k < j) already computed, and 
in this way utilize the most accurate information available. This pro- 
cedure is, in fact, easier to implement in a computer program and, 
in addition to requiring less storage, often has better convergence 
properties than Richardson’s method. It may be shown that the matrix 
representation, analogous to (46), for the ‘‘Gauss-Seidel”’ method is 
x= Wer + fe (47) 
We =F -D—-L)'u (48) 
fo = (9 —-D-L)"F, 
D, £, and U having been defined previously. 

The eigenvalue condition given in (40) is, of course, a very restrictive 
one, so that iteration procedures cannot be applied with success to 
every system of equations. However, when a convergent matrix 9% 
can be found, the methods which have been discussed offer several 
distinct advantages over a straightforward matrix inversion. For ex- 
ample, if the order of the system is N, then it can be shown that each 
iteration requires approximately N° multiply-add operations.* On the 
other hand, an inversion requires at least N° multiply-adds, so that 
the relative saving Is the ratio of the order N to the number of itera- 
tions required. It is often the case that the maximum eigenvalue is so 
small that the number of iterations required for an accuracy equivalent 
to that obtained by inversion is considerably less than N. 

Iterative methods also have the property that the solution accuracy 
is “‘adjustable’”’, in the sense that once the solution has converged to 
the point where some norm, ¢.g., €g or €, defined previously, is less 
than a specified tolerance, the iteration process can be terminated. 
This property is especially attractive in view of the fact that truncation 
errors have already been introduced, and it would therefore be super- 
fluous to accurately invert a system which is itself approximate. By 
having the option of terminating the iterative procedure, we introduce 
an additional degree of freedom by which we can optimize the computa- 
tional program. 


V. APPLICATIONS 


The iterative techniques discussed in the previous section will now 
be applied to two problems of interest, namely the H-plane discontinuity 


* A multiply-add consists of the multiplication of two complex numbers and the 
adding of the result to a third complex number. For repetitive computational 
evens, the number of multiply-adds is a measure of the computational effort 
required. 
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problem mentioned in Section IV and the analysis of TE,, — TMi: 
mode conversion at a step discontinuity in a circular waveguide. In 
both of these examples the required matrix elements may be expressed 
in convenient closed forms, which considerably reduces the required 
computational effort. 


5.1 H-plane Discontinuity in a Rectangular Guide 


In Section IV, the H-plane discontinuity of Fig. 3 was analyzed 
using an iterative algorithm which was observed to exhibit an oscillatory 
instability when initiated with a poor approximation to the actual 
solution. It was concluded that this was due to the eigenvalues of the 
iteration matrix SW being close to, or perhaps greater than unity. We 
now consider the same problem using the Gauss-Seidel method which, 
on the basis of the previous discussion, 1s expected to improve matters 
substantially. 

We assume a normalized TE,, mode incident from the smaller guide. 
Because of the symmetry of the junction, such a wave excites only 
Tk modes in both the + and — regions. The problem is, of course, 
to determine the corresponding modal coefficients for the fields on 
each side of the discontinuity. We shall present results only for the 
transmitted TE,, mode, which is the only mode propagating in the 
larger waveguide for the present dimensions, ka, = 3.5, ka, = 4.5. 

It can be shown” that the normalized vector wave functions are 
given by 


-E’ =e, Teak sin |e. (x + a) | 
1 


+B 1p : PT ,_.. 
Ey e, Py cEa: sin |e (2 + a) | ; 





(49) 


where a, , @, and 6 are the dimensions of the guide as shown and e, is 
a unit vector in the y direction. The corresponding magnetic vector 
wave functions can be found from (12). The respective propagation 


constants are 
2715 
ee a 2 Be. 
hy E (= | 


ae Ad = FE _ (ez. } 
: L 202 


From (21) and (49) we find that the coupling coefficients for the electric 
fields are given by 


(50) 


WAVEGUIDE DISCONTINUITY 667 


eal 7 = Qu a [= sin sin =a sin (2 -}- od q =) 
‘hoover Qe 2 2 |p me ay qa Z Qe 2 
Q> 


caer ae (Be sae 5) 


ay Qe 
TT Gv 
pr— 7g 


“F (51) 


The appropriate magnetic field coupling coefficients are easily found 
from (51) using the relation 


+, carat (Lees a Ae ™ y + 
ap : | a Nar ; : : ™ 
where Z, denotes the modal impedance, equal to 
CS A | gee el a 
4, = hy) (53) 


The results for the TE,, modal coefficient, “B,, as obtained using 
the Gauss-Seidel iteration method, are conveniently represented in 
Table I. Also given are the numerical values for the error parameters 
Ep, €g and €,. Truncation for this example was at 25 modes in each 
waveguide. 

We conclude that for most applications, two iterations would probably 
have been sufficient, corresponding to a saving of greater than 90 percent 
in actual execution time, compared to a matrix inversion. The reason 
for this extremely rapid convergence is, as expected, in the magnitude 
of the largest eigenvalue, which was found, using (48), to be less than 
0.078. 

As a means of establishing the convergence of the normal mode 
solution, we have plotted in Figs. 7 and 8, the mean square errors 
€, and €,, respectively, as a function of the number of modes taken. 


TABLE I— ReEsutts USING THE GAUSS-SEIDEL MErHoD FOR ANALYSIS 
oF THE H-PLANE DISCONTINUITY 


+By Energy 


Iteration Sa coefficient rms error rms error 
number Real Imaginary Ep Ep en 
1 0.97445 0.00951 —0.62X10- 0.5X1075 0.01445 
2 0.97747 0.00455 —0.74X107% | 0.49107 0.01355 
3 0.97747 0.00426 0.85107 | 0.491075 0.01350 
4 0.97747 0.00424 0.71X10-* | 0.491075 0.01349 
5 0.97747 0.00424 0.51107 | 0.491075 0.013849 
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NUMBER OF MODES 


Fig. 7— Mean square error ¢g, as function of the number of modes used, for 
H-plane discontinuity. 


These figures give genuine significance to the term ‘‘convergence in 
mean square’, since they indicate that the use of more terms leads 
to better agreement with the boundary conditions in the mean square 
sense. Again, the error is uniformly higher for the magnetic field than 
for the electric field, due to the singularity at the corner of the dis- 
continuity. 

It is finally of interest to determine the effect of a poor initial estimate 


O POINTS ACTUALLY 
DETERMINED 


MEAN SQUARE ERROR €y 





NUMBER OF MODES 


Fig. 8— Mean square error ey, as function of the number of modes used, for 
H-plane discontinuity. 
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TABLE II — RESULTS OF SPOILING ELEcTRIC FIELD At 4TH ITERATION 


TEno coefficient 


Iteration Se 

number Real Imaginary 
3 0.97747 0.00426 
4 1.94165 0.02572 
5 0.98052 —0.00063 
6 0.97747 0.00396 
7 0.97746 0.00423 


of the solution. We find that, unlike the simple algorithm discussed 
in the previous section, the Gauss-Seidel procedure is very stable, 
returning to the correct ‘‘steady state’’ solution within a few iterations 
after the spoiling was introduced. The results are shown in Table II, 
again for a truncation of 25 modes in each waveguide. 


5.2 Mode Conversion at a Step Discontinutty in a Circular Waveguide 


The second problem to which these techniques were applied is that 
of calculating the TE,, — TM,, mode conversion at an abrupt dis- 
continuity in a circular waveguide. Recent studies have indicated that 
this configuration 1s a very efficient transducer for use in dual mode 
conical horns.” The discontinuity is Ulustrated in Fig. 9 which shows 
the TE,, mode, incident from the smaller guide, being converted to 
a combination of TE,, and TM,, modes propagating in the larger guide. 

The normalized TE and TM vector wave functions are known"’ 
and, fortunately, it is possible to determine the appropriate coupling 
coefficients. or the elements of the matrix we find that 





Fig. 9—TEn — TMnu mode conversion at a discontinuity in a circular waveguide. 
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_— 2 he *hi*2J (c. 4 
| | ay aeons are (54) 
7 (Gat at)ney 
Quy" byzy E (4y =e | (Sy ) 
oe a oes ile ass a Np : 
Gig ar = Wy pee ee (55) 
(3 i W)-v (ys — DO — 1) Silty.) | 
~ Zu *nied,( x, 2) 
Pe ss pase ae es (56) 
Vin | tad 2(%2) V (yz — 1) 
oa 7 : = 0. (57) 


In (59) through (62) we have used the following notation: 


a — radius of smaller waveguide, 
b — radius of larger waveguide, 
XZ, — pth zero of J,(x), 
Yo — pth zero of J{(y). 


The elements of the matrix 3C can easily be found from the impedance 
relations of (11) and (12). 

One parameter which has been found to be useful in characterizing 
the mode conversion properties of the discontinuity is the conversion 
coefficient C’, defined as the ratio of the p-components of electric field 
for the two modes evaluated at the wall of the larger waveguide, i.e., 


vl M 


academe 
TE 
ii, 








C = 20 logio dB. (58) 


p=b 


This quantity was calculated for the particular discontinuity a = 1.05”, 
b = 1.4” over the frequency range 5.2 to 7.0 kHz, these parameters 
having been chosen for purposes of comparison with available experi- 
mental data. Truncation of the normal mode expansion was made 
after twenty-five TE and twenty-five TM modes in each waveguide. 
The iterative sequence was terminated when successive values of the 
modal amplitudes differed by less than 107°. It was found that typical 
values for the error criteria are ep & 107’, and eg, ey 0.015. These 
results indicate that for a given accuracy, a much lower value can be 
expected for ep, which is a function only of the lower-order modes, 
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Fig. 10—TEn — TMi conversion coefficient of a step discontinuity in a circular 
waveguide. a = 1.05’, b = 1.4”. 


than for eg and €y,, which depend on the higher-order terms as well. 

In Fig. 10 we have plotted the computed values of the conversion 
coeficient, defined in (58), as a function of frequency. Also shown, 
as discrete points, are experimental results obtained previously.*° The 
theoretical values are seen to be in very good agreement. 


VI. SUMMARY AND CONCLUSIONS 


In this paper, we have considered the solution of those matrix equa- 
tions which arise in the analysis of a class of waveguide discontinuity 
problems. In searching for criteria to estimate the accuracy of computed 
results, we have found that the uniqueness theorem itself yields a 
convenient set of error parameters which are easily implemented in 
the computational program. 

It is suggested that an iterative technique, particularly the ‘“‘Gauss- 
Seidel” or “point-single-step” method often leads to a rapidly con- 
verging solution, thus offering several advantages over the usual inver- 
tive procedures, particularly in the speed of execution. Of particular 
interest is the fact that when this method is applied to the analysis 
of TE,, — TM,, mode conversion at an abrupt discontinuity in a 
circular waveguide, it yields a rapidly convergent and accurate solution. 
This has been established not only on the basis of theoretical error 
criteria, but also by comparison with experimental results previously 
obtained. 
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